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Abstract 


This thesis presents a numerical model capable of simulating offshore wind turbines 
exposed to extreme loading conditions. External condition-based extreme responses 
are reproduced by coupling a fully nonlinear wave kinematic solver with a hydro- 
aero-elastic simulator. 

First, a two-dimensional fully nonlinear wave simulator is developed. The tran- 
sient nonlinear free surface problem is formulated assuming the potential theory and 
a higher-order boundary element method (HOBEM) is implemented to discretize 
Laplace’s equation. For temporal evolution a second-order Taylor series expansion 
is used. The code, after validation with experimental data, is successfully adopted 
to simulate overturning plunging breakers which give rise to dangerous impact loads 
when they break against wind turbine substructures. The impact force is quanti- 
fied by means of an analytical model and the total hydrodynamic action is finally 
obtained by adding the impulsive term to the drag and inertial ones. 

In the second main core of the thesis, emphasis is placed on the random nature 
of the waves. Indeed, a global simulation framework embedding the numerical wave 
simulator into a more general stochastic environment is developed. Namely, first 
a linear irregular sea is generated by the spectral approach, then, only on critical 
space-time sub-domains, the fully nonlinear solver is invoked for a more refined 
simulation. The space-time sub-domains are defined as the wind turbine near field 
(space) times the time interval in which wave impacts are expected (time). Such a 
domain decomposition approach permits systematically accounting for dangerous 
effects on the structural response (which would be totally missed by adopting linear 
or weakly nonlinear wave theories alone) without penalizing the computational effort 
normally required. 

At the end of the work the attention is moved to the consequences that the 
proposed model would have in the quantification of the structural risk. 
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Sommario 


In questa tesi viene proposto un nuovo modello numerico in grado di simulare tur- 
bine eoliche in mare esposte a condizioni ambientali estreme. Le simulazioni, ot- 
tenute per via diretta (i.e. a partire da azioni estreme stabilite a priori fissando 
il periodo di ritorno), accoppiano un modello numerico di propagazione del moto 
ondoso con un solutore idro-aero-elastico dell’intero sistema. 

Inizialmente viene sviluppato il modello numerico bidimensionale e comple- 
tamente non-lineare di propagazione dell’onda formulato assumendo un flusso a 
potenziale. L’equazione di Laplace ad esso associata viene risolta numericamente 
attraverso il metodo degli elementi al contorno di ordine elevato (HOBEM). L’evo- 
luzione temporale del moto viene effettuata implementando una serie di Taylor fino 
al secondo ordine. Il software, validato con dati sperimentali, è così in grado di 
riprodurre onde frangenti fino al rientro del getto. La forza di impatto, quantificata 
mediante un modello analitico, viene poi aggiunta alle componenti inerziale e di 
trascinamento in modo da stimare l’azione idrodinamica complessiva. 

La seconda fase del lavoro è dedicata allo sviluppo di un modello di simulazione 
globale finalizzato alla integrazione del suddetto solutore numerico in un ambiente 
del tutto stocastico. In particolare, dapprima una mare lineare irregolare viene gen- 
erato con approccio spettrale. Successivamente, solo su sotto-domini critici, i.e. 
nell’intorno della turbina e per intervalli di tempo in cui i frangimenti sono attesi, il 
simulatore non-lineare viene lanciato per un’analisi più raffinata. Questa strategia 
di decomposizione del dominio permette di mettere in conto in modo sistematico 
quegli effetti impulsivi, che verrebbero altrimenti ignorati dai modelli lineari, senza 
penalizzare lo sforzo computazionale normalmente richiesto. 

Infine, il lavoro si chiude mostrando le ripercussioni che il modello proposto può 
avere in termini di affidabilità e sicurezza strutturale. 
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Kurzfassung 


In dieser Arbeit wird ein numerisches Modell zur Simulation von Offshore-Wind- 
energieanlagen unter extremen Lasteinwirkungen entwickelt. Dazu wird ein voll- 
ständig kinematisch nichtlineares Wellenmodell mit einem hydroaeroelastischen Mod- 
ell kombiniert. 

Zunächst wird das instationäre nichtlineare Problem der freien Wasseroberfläche 
unter Verwendung der zweidimensionalen Potentialtheorie beschrieben. Die sich 
ergebende Laplace-Gleichung wird mit einer Randelementmethode höherer Ord- 
nung räumlich diskretisiert. Für die zeitliche Entwicklung wird eine Taylor Reihe 
zweiter Ordnung verwendet. Nach Abgleichung mit experimentellen Daten wird 
der entwickelte Algorithmus angewendet, um die für die Stoßbelastung von Wind- 
kraftanlagen ursächlichen überschlagenden brechenden Wellen zu simulieren. Die 
gesamte hydrodynamische Last wird schließlich durch ein analytisches Modell be- 
schrieben, bei dem ein Term, der die Stof&wirkung der Wellen berücksichtigt, zu den 
Längs- und Trägheitskräften hinzugefügt wird. 

Im zweiten Teil der Arbeit wird das Wellenmodell in eine Simulationsumgebung 
eingebettet, welche die stochastischen Natur des Wellengangs erfasst. Hierbei wird 
zuerst ein linear beschriebener breitbandiger Seegang mithilfe des Spektralansatzes 
erzeugt. Beschränkt auf kritische Bereiche in der räumlichen und zeitlichen Simu- 
lation wird im Anschluss das vollständig nichtlineare hydrodynamische Modell für 
eine genauere Lösung herangezogen. Die kritischen Bereiche sind auf die nähere 
Umgebung der Windenergieanlage beim Eintreffen der brechenden Welle begrenzt. 
Diese Substrukturtechnik erlaubt es, fiir die Strukturantwort mafgebende Effekte 
systematisch zu erfassen, die bei einer Verwendung von linearen oder schwach nicht- 
linearen Wellentheorien komplett vernachlässigt werden, ohne dabei den herkömm- 
lichen Rechenaufwand substantiell zu erhöhen. 

Zum Abschluss der Arbeit wird diskutiert, wie sich das vorgestellte Modell auf 
die Quantifizierung des Risikos der Struktur auswirkt. 
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Chapter1 
Introduction 


This introductive chapter presents the framework of the thesis and outlines the moti- 
vations and the organization of the work. 

At the beginning it also tries to give a global view about the wind energy market with 
the aim of pointing out that research on (offshore) wind energy is nowadays truly crucial. 


1.1 Wind energy 


To introduce the topic of this thesis and to get a broader idea about the global 
energy system, it could be useful to start with some preliminary notions about the 
current worldwide and European energy consumption status as well as the objectives 
in matter of greenhouse gas emissions!. 

Only starting by clarifying the objectives and, above all, the departing point 
for each country, it is possible to understand how important is nowadays to invest 


research sources in wind energy. 
1.1.1 European and world energy scenario 


During the last decade with no doubt wind energy has been representing the 
leading renewable energy source. And, according to government plans, it will keep 
being the leading renewable source for many years. 

The worldwide energy scenario is represented in figures 1.1 and 1.2 which show 
the global annual wind power installed capacity in the period 1996-2008 and the 
annual wind power capacity installed by region, respectively. 

Wind power is the fastest growing power generation technology in the EU with 
more than 35% of all new energy installations in 2008. It is also interesting to note 
from figure 1.2 that European wind energy installation has been leading the global 
installation since 2003. 

Figure 1.3 shows that in 2008 wind energy installation was definitely dominating 
other energy sources. Indeed, only in this year, upon a total installed capacity of 
23.851 GW in Europe, approximately one third was represented by wind anergy. 

Almost 8.9GW of new wind turbines installed in 2008 brought European wind 
power generation capacity up to nearly 66GW. Another promising sign, see fig- 
ure 1.4, is the diversification of the European market. 2008, in fact, saw a much 
more balanced expansion with not negligible contributions given by Italy, France 
and the UK. 


1Data and statistics here presented are all update to 2008. Currently, official data from EWEA 
and GWEC for 2010 are not yet available. 
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GLOBAL ANNUAL INSTALLED CAPACITY 1996-2008 


1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 
1,280 1,520 2,520 3,440 3,760 6,500 7,270 8,22 8,207 11,521 15,245 19,865 27,051 


Figure 1.1: Global annual wind power installed capacity, 1996-2008. 


ANNUAL INSTALLED CAPACITY BY REGION 2003-2008 
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Figure 1.2: Annual wind power installed capacity by region, 2003-2008. 
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* Geothermal, peat and waste 
** This Is a preliminary figure for solar photovoltaic installations (source: European Photovoltaic Industry Association (EPIA). 


Figure 1.3: New power capacity installed in Europe in 2008. 
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However, Germany continues to be Europe’s leading market, both in terms of 
new and total installed capacity. To confirm this, figure 1.4 shows that over 1.6 GW 


of new capacity was installed in 2008 and this brought the total German capacity 
up to nearly 24GW. 


TOP 10 TOTAL INSTALLED CAPACITY 2008 TOP 10 NEW CAPACITY 2008 
Rest of world Rest of wortd 
us us 
Portugal P Mrs ay 
Denmark & & 
UK 
ps France 
France fim fe 
naty italy 


P |- -> 
india 


Spain India 
Mw % Mw % 
USA 25,170 208 USA 8,358 30.9 
Germany 23,903 19,8 China 6,300 233 
Spain 16754 139 India 1,800 67 
China 12,210 10.1 Germany 1,665 62 
India 9,645 80 Spain 1,609 5.9 
Italy 3736 31 Italy 1,010 37 
France 3,404 28 France 950 35 
UK 3,241 27 UK 836 31 
Denmark 3,180 26 Portugal 712 26 
Portugal 2862 24 Canada 526 19 
Rest of world 16493 13.8 Rest of wodd 3,285 122 
Total top 10 104,104 862 Total top 10 23766 87.8 
World total 120798 100.0 World total 27,051 100.0 


Figure 1.4: Top 10 global capacity installed, total and in 2008. 


As reported in [15], it is also worth mentioning that among the growing European 
markets in 2008, Italy experienced a significant leap: over 1 GW of new wind turbines 
came on line in 2008, bringing the total installed capacity up to 3.7 GW. 


A more detailed map about the cumulative installed capacity state by state is 
given in figure 1.5. 


| Installed in 2008 | Cumulative, end of 2008 
Total EU-27 | 8484 MW | 64935 MW 


of which Offshore | 357 MW | 1471 MW 


Table 1.1: European wind power capacity. 


As shown in table 1.1, at the end of 2008, there were 65 GW of wind power ca- 
pacity installed in the EU-27 producing 142 TWh hours of electricity which satisfies 
4.2% of the whole EU electricity demand. This means that at the moment offshore 
wind energy is able to satisfy only 0.1% of the whole EU demand. This datum 


makes more understandable how challenging are the targets fixed by EU which will 
be shortly recalled in the next section. 
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Wind power installed in Europe 
by end of 2008 (cumulative) 


Faroetislands 
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European Union: 64,935 MW 
Candidate Countries: 452 MW 
EFTA: 442 MW 

Total Europe: 65,933 MW 


Figure 1.5: State by state cumulative installed capacity at the end of 2008. 


In 2008 US wind industry was able to install 8.36 GW marking an increase in 
generating capacity of 50% in a single calendar year. The 2008 US growth repre- 
sented about 42% of new electricity generating capacity added in the United States 
during the year, establishing wind as a mainstream energy source for the country 
(second only to natural gas) in new generating capacity. US total wind generating 
capacity in 2008 was more than 25.17 GW, producing enough electricity to power 
the equivalent of close to 7 million households and to meet over 1% of total US 
electricity demand. 


1.1.2 Short and long term objectives 


Focusing on the European situation, in the Strategic Research Agenda (SRA) 
- a document prepared by the Wind Technological Platform (TPWind) in 2008 - 
fundamental objectives in matter of wind energy development have been fixed. They 
are divided into: 


Short term targets: within 2020 reduction of greenhouse gas emission by 20% 
and ensure 20% of renewable energy sources in the EU; 


Long term targets: decarbonization, 60 - 80% reduction of the greenhouse gas 
emission. 
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To meet the 2020 targets, among many other research lines, for the European 
Commission it is imperative to: 


“Double the power generation capacity of the largest wind turbines, with 
offshore wind as the lead application.” 


In particular for the offshore wind, the Strategic Research Agenda establishes 
the following special objectives to be achieved within 2030: 


e More than 10% of Europe’s electricity should come from offshore wind; 


e Make the offshore generating costs competitive with other sources of power 
generation; 


e Make commercially mature the technology for sites at any distance from shore 
with a water depth up to 50m: 


e Full-scale proven technology to dominate deep-water sites. 


Moreover, together with the above targets, five research topics have been prior- 
itized: 


e Substructures; 

e Assembly, installation and decommissioning; 
e Electrical infrastructure; 

e Turbines; 

e Operations and maintenance; 


With respect to on land standard designs, the offshore environment does intro- 
duce significant additional elements which have to be carefully considered, especially 
in designing the support structures. Knowledge about modeling the wind and ro- 
tor aerodynamics developed for onshore sites are generally enough and do not need 
deep changes when moving in the offshore environment. Some adjustments are made 
just due to the different wind characteristics of offshore sites (e.g. strong difference 
in the roughness length and turbulence intensity). Figure 1.6 gives an example of 
different wind shears for on- and off-shore sites, respectively. 

On the contrary, for offshore plants, the concept of support structure has to be 
entirely rearranged. For this reason research on the substructure is always prior- 
itized both directly, by improving the technology itself, and indirectly, that is by 
developing more accurate models to estimate the combined wind-waves action. In 
fact, in addition to the previous research topics, the Strategic Research Agenda 
establishes also the following priorities: 


e Development of fully integrated wind-wave-current interaction models; 
e Development of new substructure concepts; 


e Development of improved design methodologies to extend the life of structures, 
to reduce costs and to incorporate risk-based life-cycle approaches. 
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Figure 1.6: Typical wind shears for land and offshore sites. Figure from [1]. 


Consider that only the substructure represents approximately 25% of the whole 
investment and if we think that forces (and consequently costs) used to design the 
substructure increase with the square of wind/water velocity, then it appears clear 
the importance of accurate models for loads prediction. 

The objectives described above for wind energy development are based on the 
central fact that Europe has a remarkable wind potential. Figure 1.7 shows the map 
of onshore wind potential. Considering that the minimum value of the mean wind 
speed to make cost-effective a wind power plant is approximately 4m/s, it results 
that most of the European areas possess a wind energy potential. 

The offshore potential is depicted in figure 1.8, which shows that in addition 
to the North and Baltic Seas, also some Mediterranean areas, for example between 
Greek and Italian coasts, the wind resource can be exploited. 
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Figure 1.7: Onshore wind potential. European Wind Atlas. Copyright 1989 by Risg- 
National Laboratory, Roskilde, Denmark. 
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Figure 1.8: Offshore wind potential. European Wind Atlas. Copyright 1989 by Risg- 
National Laboratory, Roskilde, Denmark. 
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1.2 General nomenclature 


Before entering further into the topic of this thesis, it would be useful to pro- 
vide some general terminology which will be frequently used throughout the text. 
Figures 1.9 and 1.10 help in this task. 


Rotor shaft 
and bearings 


Gearbox 


Rotor hub and Rotor brake 


blade pitch mechanism | Generator 
ee h 


a & |} Elektrical switch bo 
== and control systems 


Rotorblade Bedplate 


Yaw system 


Power cables 


—— Tower 


Grid connection (transformer) 


Foundation | U 


Figure 1.9: Main components of an horizontal axis wind turbine. Figure from [1]. 


Figure 1.9 in particular shows the essential components of the upper part of the 
wind turbine. It is referred to an onshore case with a superficial foundation. On the 
contrary, figure 1.10 provides more details about the substructure, which is defined 
as the structural subpart included between the sea bed and the platform. Among 
the three types of substructures sketched in the figure, the so called “monopile”, the 
first from the left, is the one supporting the baseline reference model adopted in 
this work. 
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Figure 1.10: Main components of the support structures of an horizontal axis off- 
shore wind turbine. Figure from IEC61400-3 [2]. 


Figure 1.11: Offshore wind farm Utgrunden off the southern Swedish North-Sea 
coast (7 wind turbines of 1.5 MW each). 


An example of a wind farm made of monopile supported wind turbines is shown 
in figure 1.11. 
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1.3 Modeling offshore wind turbines 


Together with the general motivations regarding the global need in boosting 
wind energy production, there are some more technical lines of reasoning which 
justify the research on offshore wind turbine. 

Offshore wind turbines are very sophisticated systems that can be only analyzed 
by adopting integrated multi-physic models. There are four coupled disciplines in- 
volved: aerodynamics, structural dynamics, hydrodynamics, soil dynamics. Each of 
them relates to the rotor, the tower, the substructure, the foundation, respectively. 
Figure 1.12 tries to sketch this concept, indeed, moving from right to left, the four 
isolated subjects are applied to the four main parts of an offshore wind turbine 
and thus coupled into a unique system which should render the reality as much as 
possible. 


I 
I 
' 
I 
I Eructurai dynamics 
I 
' 


~i . 
Substructure 3 Hydrodynamics 


Soil dynamics 


Figure 1.12: Coupled disciplines in a unique system. 


The current standard technique to analyze and design offshore wind turbines 
starts from a real structure (or from a scheme if one is to design it) and, once all 
the structural and mechanical properties are known, the designer collects the system 
and environmental variables so that it is possible to provide input data for adequate 
numerical simulations. 

In the present case, the time domain solver used is FAST, which gets the aero- 
dynamic loads acting on the rotor blades by invoking AeroDyn, see scheme in fig- 
ure 1.13. Both solvers have been developed at National Renewable Energy Labora- 
tory (NREL). In some cases, however, solvers for simulating wind turbines imple- 
ment some additional hydrodynamic routines which provide the wave kinematics 
and accordingly the hydrodynamic forces to reproduce also the offshore environ- 
ment. The linear wave theory, together with Morison’s equation, is commonly im- 
plemented and this permits to simulate most of the condition that a wind turbine 
may experience during operation. In fact, looking at the green rectangles in fig- 
ure 1.13, it is clear that when the aim is to predict state of failure induced by long 
term fatigue accumulation, or more in general for all those cases characterized by 
relatively small wind speeds, this procedure is sufficient. 

On the contrary, there are design load conditions characterized by extreme wind 
which in turn generates severe highly nonlinear waves that cannot be modeled by the 
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Figure 1.13: Commonly adopted scheme for offshore wind turbines simulations. 
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wind-waves actions. 
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linear theory anymore. Later on it will be also shown that for really extreme wind 
conditions also the most common (weakly) nonlinear wave theories are not sufficient 
when phenomenal events, such as plunging breaking waves, occur. Therefore, the 
simulation scheme alone sketched in figure 1.13 is no longer valid. 

For this reason the present thesis attempts to give a contribution to fill the 
gap, so that, as highlighted by the red rectangles of figure 1.14, an improved hy- 
drodynamic model is provided with the aim of setting up and advanced numerical 
simulator able to reproduce not only irregular linear random seas associate with 
the normal wind conditions, but also extreme winds which can give rise to very 
dangerous and destructive phenomena for an offshore wind turbine such as impacts 
due to plunging breakers. 


1.4 Structure and scope of the thesis 


The thesis is articulated in four main parts. After this introductory chapter, 
which attempts to point out both the general and scientific motivations of the 
whole work, firstly the general Risk Management Chain is tailored on the specific 
features characterizing offshore wind systems. The crucial point which is stressed 
in this context lies in the fact that any design of wind energy converters must 
always assure a certain minimum system reliability level. Therefore, to estimate the 
probability of failure of such systems, a systematic procedure needs to be applied. 

Chapters 3, 4, and 5, as sketched in scheme 1.15, are devoted to the development 
of a new numerical code whose importance in terms of structural safety will be 
remarked in the concluding Chapter 6. 


What is the Global wind energy NEO 
scenario? Why research in this field? da i 
Why to fit the general Risk Man- Risk Manage- 
agement Chain to offshore wind «> ment Chain Aeroelastic Solver 
turbines? What are the needs? (Ch. 2) 


Is it possible to improve numer- s 
ical models without penalizing Re rical Core > BEN, for ul 
the computational costs? How? ERD Nonlinear Waves 


Which are the benefits of the pro- mese 
posed numerical tool in assessing the $ 
pae g effects on Risk 
probability of failure of the system?<— 
And on the general quality of OWTs 
model? Are the needs stated above met? 


Coupled Time 
Domain Solver 


Assessment 
(Ch. 6) 


Figure 1.15: Scheme of the work. Blue: state of the art, red: development of a new 
numerical tool; green: initial targets and improvements evaluation. 


Figure 1.15 helps to understand the structure of the work. The central column 
denotes the main parts of the thesis. Each of these parts tries to provide an answer 
to the key questions stated on the left-hand side. The blue elements in scheme 1.15 
denote either what concerns the state of the art or anything which has been only 
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used in the present work but not investigated. The red parts (the ellipse and the two 
related sub-cores), on the contrary, refer to parts which have been either entirely 
developed (as the case of the Boundary Element code for fully nonlinear water waves 
simulation) or improved and embedded into a new design tool. Finally, the green 
rectangles represent the motivation and objectives achieved. In particular, Chap- 
ter 2 rises some issues concerning the structural risk of offshore wind turbines, then 
Chapter 6, by using results of the central core, attempts to answer the questions. 


Chapter2 
The risk management chain of offshore wind turbines 


This chapter has the scope of applying the general probabilistic risk management chain 
to investigate the probability of failure of offshore wind turbines. 


2.1 Cost and structural safety 


In accordance with the new research lines in the field of offshore wind energy 
emphasized in the introductory chapter, fully probabilistic based design procedures 
are more and more necessary and desirable. However, it must be remarked that even 
the most sophisticated probabilistic model will return unsatisfactory results when 
the physical model being simulated is not well reflected in the numerical idealization. 

As already anticipated, the major needs for offshore systems today are repre- 
sented by novel design procedures, e.g. structural optimization, which assure the 
minimization of the costs under the constraint, among others, of a fixed Structural 
Reliability (SR). The latter is regarded as the probability that the structure under 
consideration exhibits a proper performance throughout its lifetime. The defini- 
tion itself intrinsically contains a very sensitive point: what does it mean “proper 
performance” A closer look at this issue will be taken in the next section. 

Minimizing costs under a fixed Structural Reliability level can be thought as a 
constrained optimization problem 


Cost (X) + min (2.1) 
SR (X) > SR* (2.2) 


where X is the vector collecting all design variables involved both on the loading side 
and the structural strength side, while S R* is the minimum acceptable structural 
safety level. 

Establishing the lower structural reliability level SR* mostly concerns cost- 
benefit analysis and lies beyond the goals of this work. On the contrary, the opti- 
mization problem shown in equations (2.1) and (2.2) highlights how strongly the 
cost reduction is linked to a proper estimation of SR. This concept results even 
clearer by recalling the general principle that the better the SR is evaluated, the 
larger its expected value [16]. 

Therefore, since SR is never an absolute measure - it considerably depends 
on many factors, such as the accuracy level of the idealized structural model, the 
number of uncertainties and their statistics, etc. - in this thesis the attention is 
most paid to restrict as much as possible the sources of uncertainties solely to those 
parameters which are intrinsically random (aleatory uncertainties) by providing a 
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more accurate numerical tool to account for the extreme hydrodynamic actions 
due to overturning breaking waves. It will be presented, in fact, a numerical model 
which, without assuming any a priori restriction about the magnitude of water 
waves nonlinearity, reduces the model uncertainties, especially those deriving from 
modeling the loads. 


2.2 The risk management framework 


More analytically, it is possible to quantify the structural reliability as SR = 
1 — Py where Py denotes the probability of failure of the system, or the Element at 
Risk (EaR), under consideration. Hence, the issue opened above concerning when it 
is possible to say that a structure is exhibiting a proper performance and when not, 
has now simply been shifted to determine the probability of failure of the system. 

The difficult task of determining the probability of failure, especially when the 
system is very articulated as in the case of offshore wind turbines, can only be 
carried out by invoking a proved and systematic approach which is usually referred 
to as Risk Management Process [3], [17], [18], [19] [20]. 

As already mentioned, such a general approach needs to be adapted depend- 
ing on the specific application. In fact, it might happen that not all the phases 
characterizing the Process are of primary interest for the current application. 

The first level subdivision of the whole risk management framework is made up 
of the three following crucial steps [3]: 


1. Risk Identification; 
2. Risk Assessment; 
3. Risk Treatment; 


The first step, Risk Identification, is decisive for all the remaining procedure. 
In this phase the system must be declared and the sources of all possible hazards 
menacing the system should be carefully identified. In the present case the system 
could be represented by a group of components of which an offshore wind turbine 
is made. For example, one could focus on the rotor blades, the nacelle (including 
all the electrical and mechanical components such as generator, gear box, brake 
system, etc.), the tower and finally the foundation. 

However, in the case being investigated in this thesis the system is represented 
by the tower and the performance we are interested in is its capability of safely 
carry all the loads acting on the entire structure. This choice is justified by the fact 
that since extreme environmental conditions are simulated, the rotor is always in 
the parked condition so that an ultimate load case condition will always involve 
primarily the support structure!. 

The threat which may endanger the system is represented by a combined wind- 
waves action. As in this work we only consider wind-generated waves, in the hazard 
identification we distinguish a “driver” (main) hazard, represented by the wind and 
a “driven” (induced) hazard represented by the waves. Note that this distinction 


In future the effects of the foundation should also be considered, but in the current application 
the monopile is considered rigidly connected to the sea bed. 
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Figure 2.1: Schematic representation of the two main loading conditions an offshore 
wind turbine may experience. Different failure types have to be investigated with 
different tools. 


by no means should lead to believe that waves are less dangerous or secondary in 
terms of load intensity compared with the wind. 

These two hazards are not statistically independent and thus in the subsequent 
phase, Risk Assessment, the correlated multi-hazard scenario needs to be carefully 
analyzed. 

The Risk Identification phase ends only when it is possible to answer the fol- 
lowing question: given the system, what is the possible state of failure that it may 
experience? Figure 2.1 helps to answer this question. Indeed, the first most impor- 
tant distinction has to be made between ultimate failure condition and long-term 
damage accumulation failure. What is crucial to stress is that not only different 
statistics are involved depending on the failure type under investigation, but also 
the way of modeling the actions has to be be properly adjusted. 

It is known that offshore structures are basically exposed to wind, wave and cur- 
rent loads, see the schematic representation in figure 2.2, and they may fail because 
of two different loading conditions: (i) extreme wind coupled with highly nonlinear 
extreme waves, which lead to the so called ultimate failure condition and (ii) long 
term fatigue accumulation that gives raise to the growth of cracks |21]. In general 
the fatigue failure results from the accumulation of damage induced by fluctuat- 
ing loads. When the material is exposed to a continually changing internal state 
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Figure 2.2: Schematic representation of the main loading actions on an offshore 
wind turbine. 


of stress, it slowly deteriorates initiating cracks which will eventually lead to the 
material rupture. Such failures caused by damage accumulation are not investigated 
in this thesis. Contrarily, we will only focus on failures induced by extreme loading 
conditions. In particular, it will be assumed that our system reaches a state of fail- 
ure whenever the state of internal forces due to the combined action of wind and 
waves attains a certain fixed value above which the bearing loads capability of the 
tower is compromised?. 

A deep understanding of the failure scenario paves the way for the next phase 
of the risk management framework: the Risk Assessment. 


Risk Assessment 


The Risk Assessment step represents the core of the whole risk management 
process and it is divided in two subparts: 


e Risk Analysis; 
e Risk Evaluation; 


The above approach is referred to the model proposed in [3], see figure 2.3. 

This phase aims at quantifying the risk. To this end it is of primary importance to 
employ as much accurate as possible prediction models to estimate both the hazard 
intensity and its frequency of occurrence. The risk is quantified by the following 
expression 


Risk = Probability of failure x Losses [Losses unit /time] 


? Additional loads might be provided by ice impacts and earthquakes, but in this context they 
are not taken into account. 
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Figure 2.3: The Risk Assessment phase. Image from [3]. 


where Losses indicates all the expected consequences which may happen due to 
the failure of the system. Typical losses are: economical, fatalities, effects on the 
environment, etc. The estimation of consequences given the failure of the system 
remarkably relates to the nature of the system. Consequences associated with a 
failure of support structures of offshore wind turbines may be very significant in 
terms of monetary loss, while fatalities have less probability of occurrence than in 
the case of structural failures of civil buildings. Also the impact on the environment 
would not be dramatic as a failure of nuclear power plants. For these reasons, in 
the present work the risk analysis is restricted only to the first term of the risk, 
that is the the probability of failure of the system. The consequences - and thus the 
quantification of the losses in case of failure of the system - is not part of this thesis. 

The probability of failure involves the Vulnerability of the system, the Hazard 
and the Exposure. According to [22], [23], [24], and assuming that the uncertainties 
associated with the loads are definitely dominating those related to the structure 
[25], for the specific case under consideration, the probability of failure can be 
expressed as follows 


Pig = P (LS) = / P (LS|IM) p (IM) dIM (2.3) 


where LS denotes a limit state which in this case it is assumed to represent both 
the measure of the damage induced by the hazard and the structural response, 
while IM is the intensity measure of the hazards. In the light of equation (2.3) it 
is possible to understand the following conceptual equation for the probability of 
failure 


Prait = Vulnerability x Hazard x Exposure 


where the Vulnerability is meant as the probability of a certain structural damage 
(or structural response) given the intensity measure of the hazard: P (LS|IM), while 
the Hazard is meant as the probability distribution the hazard intensities: p (IM), 
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finally the Exposure is intended as the probability that such a Hazard meets the 
system. 

As already mentioned, in this work we are mainly interested in the bearing load 
capacity of the substructure so that the structural response is quantified by means 
of the tower base bending moment M,; and, as a consequence, the Limit State 
function which identifies the state of failure is given by 


Limit State (LS): My: > My, & Failure 


where M7, is the ultimate resistant tower—base bending moment. 

It has to be noted in equation (2.3) that the hazard intensity measure is a 
vector (typed in bold). Uncertainties from the load side are provided by several 
variables: the mean wind speed, the turbulence intensity, the wind shear exponent, 
yaw misalignment, significant wave height, zero up-crossing wave period, current 
speed, etc. However, in agreement with what observed in [25], the randomness in 
the load model is restricted only to the three variables which have the greatest 
influence on the response of support structures. These are: the mean wind speed U, 
the significative wave height H, and the zero up-crossing wave period T,. Thus, by 
setting IM = (U, H,, T), equation (2.3) becomes 


Prait=P (My: > Mj,) = fP (My > M}.|U, Hs, T;) p (U, Hs, T.) AUAH,dT, 


(2.4) 
where £ is the environmental domain which represents all the possible variations of 
U, H, and T,. 

The conditional distribution of the response given the set of sea state parameters 
is usually referred to as short-term distribution. By multiplying the conditional 
response times the joint probability density function of the sea state parameters U, 
H, and T, and integrating over the domain £, the so called long-term distribution 
of the response is obtained. 

According to [25] and as sketched in figure 2.4, the procedure to get the prob- 
ability of failure conditioned on the environmental parameters is articulated in the 
following main steps: 


e Definition of a joint wind-waves probabilistic model (from measurements or 
hindcast data); 


e Generation of wind and waves time histories from which aerodynamic and 
hydrodynamic loads are derived; 


e Time domain simulation and analysis of response time series; 


At this point some statistics of the time series permits to obtain the distribution 
of the response conditioned on a given set of environmental parameters. To ob- 
tain such a conditional distribution several methods can be applied, for a detailed 
description we refer to [25], [26]. 

As already mentioned the long-term distribution is then obtained by the integra- 
tion over the domain of all sea state. Finally, to extrapolate the response distribution 


3Here this value is assumed to be deterministic. 
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Figure 2.4: Schematic representation to obtain the short-term response given the 
environmental parameters intensity. 


Base moment 


for longer period, e.g. 50 years, due to the hypothesis of independent sea states (once 
the sea state duration has been established), it is possible to calculate the number 
of independent sea state in a fixed time period and thus the probability distribution 
of the response for such a return period. 


Finally, the probability of failure, i.e. the probability that in such a period of 
time the maximum response is higher than the ultimate strength, can be computed. 


The methodology above described is known as “response—based” approach. Al- 
though it represents a rational method to estimate the probability of failure of the 
system, it has some disadvantages: it requires a considerable computational effort, 
it cannot be used in the early stage of the design; it usually adopts the same load- 
ing model independently on the intensity of the environmental parameters. Thus, 
no possibility of accounting the increasing waves nonlinearity as the wind speed 
increases. 


An alternative procedure used to estimate the extreme response of offshore wind 
turbines is the so called “external conditions—based” approach. In this case it is 
assumed that the extreme response is induced by the extreme loads. In this thesis 
this methodology is used because our scope is to capture the response when the 
structure system is exposed to extreme events such as breaking waves. It will be 
shown in fact that, according to IEC61400-3, when 50-year return period storm 
are simulated, then breaking waves occur causing dangerous impacts against the 
substructure. 


Within this second approach, first the joint probability distribution for the three 
main random variables U, H, and T, has to be built [21], [11], [27]. Then from the 
joint model the environmental contour may be obtained by fixing the requested 
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return period. Next, assuming that the turbine is parked or standstill!, and that 
the extreme response is always proportional to the wind speed and the significative 
wave height, the triplet of environmental parameters which induce the maximum 
response can be easily found among the infinity triplets of U, H, and T, laying on 
the environmental contour [28], [21], [27], [26], [29], [30]. 

For the practical application in this thesis, however, because of lack of measured 
data and also because it is not central for our purpose, the extreme environment, 
is established in a simplified way, that is by first setting the wind intensity in 
agreement with the extreme turbulent wind speed model (EWM) as in IEC61400-1, 
then by deriving the sea severity - associated with such a wind model - by using 
values for H, and T, recommended by literature [10], [31], [32], [11], [33], among 
others. Details on the wind—waves loading model are discussed in section 5.4. 

The point we want to to stress in conclusion of this chapter is that the loading 
model employed to predict the system response needs to be adjusted according to 
the intensity measures of the hazard. Currently, appropriate numerical tools capable 
of systematically predicting the structural response when the hazard is characterized 
by fully nonlinear waves seem to lack in literature. 


4In particular we will only focus on the substructure bearing load capacity rather than on the 
internal forces of the rotor blades. According to what said in the previous section, this restriction 
is reasonably justified because in case of extreme scenarios it will be assumed that the turbine is 
parked, so that blades are set with 90° pitch angle in order to be not affected by any lift force. On 
the contrary, in extreme environmental conditions tower and substructure result mostly involved. 


Chapter3 
Aerodynamic model 


In this chapter the basic notions concerning the rotor blade aerodynamics on which the 
solver used in the next chapters is based are shortly recalled. The key concepts of the Blade 
Element Momentum are presented step by step starting from the disc theory. 


3.1 Basics on wind turbines aerodynamics 


The primary goal of a wind turbine is to subtract kinetic energy from the wind 
to transform it first into mechanical energy and then into electrical energy. The 
conversion of the wind kinetic energy into mechanical energy takes place when the 
air flows through the rotor disc. Given the upstream air flow velocity and making 
use of some basic fluid dynamics, it is possible to calculate both the velocity at 
the rotor disc and in the wake, provided that the axial flow induction factor a is 
known. The maximum achievable value of power coefficient is known as Betz limit 
and represents only a theoretical value. What so far introduced is usually named 
Momentum Theory and it is marked by the fact that it does not employ neither any 
event which occurs locally at rotor blades nor the shape and the number of blades. 

To calculate the torque and power developed by the rotor a more sophisticated 
model involving lift and drag forces on the blades is adopted. To this aim, first 
the axial wind velocity at the disc is composed with the tangential velocity, which 
depends on the rotor angular velocity as well as on the tangential flow induction 
factor a’, then, given the aerodynamic coefficients of each elemental segment of the 
blades, it is possible to compute the drag and lift forces. An iterative procedure 
permits to calculate the induction factors a and a’ which finally lead to know the 
torque and, as a consequence, the power developed by the rotor. The latter, divided 
by the maximum available power, gives the expression of the power coefficient. Cp. 
The variation of Cp versus the tip speed ratio represents the performance curve of 
a wind turbine. 

In the following section details about the method will be presented and the 
meaning of the terminology here introduced will result clear. 


3.2 Momentum theory 


Let us start from the stream tube concept sketched in figure 3.1. The mass of 
fluid m passing a generic cross section A of the stream tube is given by m = pV, 
where p is the air density and V the volume of fluid. The mass flow rate, that is the 
volume of fluid passing the cross-section A per unit time is just given by Q = AU, 
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where U is the velocity of the air particles passing through the considered section!. 
3.2.1 Axial momentum 
Of course, mass flow conservation imposes that 
Quo = Qa = Qu (3.1) 


where the subscripts oo, d and w - according to figure 3.1, denote the far upstream 
undisturbed wind velocity, the disc, and the wake wind velocity, respectively. 
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Figure 3.1: Stream tube. Wind Energy Handbook, Burton et al. Wiley 2001 [4]. 


Equation (3.1) can also be written as follows 
AU = Ada = Aww (3.2) 


Newton’s second law can be naturally applied at the disc by introducing the 
momentum K = mU. At the disc, the momentum rate writes as Kg = pAaUaU. 
Thus, throughout the disc the momentum rate experiences an overall change given 
by 

AKa = pAqUaAU = pAqUa (Uso — Uw) (3.3) 


and consequently the total forces acting at the disc, basically given by the pressure, 
must equal 


Aa (pt — py) = pAaUa (Uso — Uw) (3.4) 


Now, by introducing the azial flow induction factor a, which permits to express 
the velocity at the disc through the far upstream undisturbed velocity U» as Ug = 
Us (1 — a), equation (3.4) turns into 


Apa = pU» (1 — a) (Uso — Uw) (3.5) 


1Q = £ (Ads) = Ad = AU, where ds is the disc thickness: the distance the air particles would 
cover in a time interval of dt. 
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The pressure gradient Apg can be computed by using twice Bernoulli’s equation 
(once between the upstream undisturbed section and the disc section and once 
between the disc and the wake sections). This leads to write 


1 
Apa = 5p (U2, - U2) (3.6) 
that replaced into equation (3.4) finally gives 
Uw = (1 — 2a) U> (3.7) 


Next, by making use of equation (3.7) into (3.4), the total force acting at the 
rotor disc (the thrust) is given by 


Fy = 2pAqU2,a (1 — a) (3.8) 
and, accordingly, the power developed is 
Prieta = FiUa = 2pAqU3,a (1 — a)? (3.9) 


It is straightforward now to see that from the all available power in the wind, 
machines can only extract a share given by the so called power coefficient Cp defined 
as follows 


= = 4a (1 — 3.10 
= Tea Tela (3.10) 

The above coefficient is maximum when a = 1/3, therefore we have 
Chmax = 9.593 (3.11) 


that is known as Bet’s limit. It represents an ideal value which proves that the most 
exploitable energy is about 60 % of the available wind power. 


3.2.2 Angular momentum 


We are still not considering the blades aerodynamics. Since the final scope is to 
get the aerodynamic forces acting on the blades, it is necessary to go more in depth 
by adding some concepts about the angular momentum. 

We assume that an air particle past the rotor disc has a tangential velocity 
U; = 2a'OQR where Q is the angular velocity of the rotor, while a’ is called tangential 
flow induction factor. Refer to figure 3.2. 

It is easy thus to derive that, given an elemental volume of air, the infinitesimal 
angular moment writes as $K "9 = fF x dmU; = dm2a'Qr, where the elemental mass 
is ôm = provor. By integrating along the circumference and taking the derivative 
with respect to time, we get 


6K "I = pS AqUoo (1 — a) 2a' Dr? (3.12) 


which represents the change of angular momentum at the rotor disc regarding an 
elemental annulus. Therefor, it is straightforward to obtain an alternative expression 
for the power yield? 


bPyieta = PO AaUco (1 — a) 2a'Q? r? (3.13) 


?Power is given by the couple times the angular velocity. 
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Equating equation (3.9) with the above result we get 
a’? = a(1—a) (3.14) 


where A, = Qr/Us. is the local speed ratio. 

Throughout the rotor disc thickness the tangential velocity of an air particle 
varies from zero (upstream) to 2a’Qr (downstream), so that in the middle a linear 
model permits to consider that the tangential velocity is a’Qr. See figure 3.2. The 
latter, along with the tangential velocity of the blades, gives rise to a net tangential 
flow experienced by the blade element equal to Qr + a'Qr = Or (1 +’), therefore 
the total relative velocity experienced by the blade is the following. See figure 3.3. 


W = uz (1-0)? + 92r? (1+ a)? (3.15) 


2a'Qr 


Pp—3p(2a'Qr)? 


Rotor motion 


Figure 3.2: Tangential velocity growing across the rotor disc thickness. Wind Energy 
Handbook, Burton et al. Wiley 2001 [4]. 


Furthermore, from figure 3.3 it can also be readily set 


W sind = Ux (1 — a) (3.16) 
W cos¢ = Or (1 — a’) (3.17) 


where 6 = a + p. The angle 8 amid the airfoil zero lift line and the plane of the 
rotor disc is named pitch angle, while a is the angle of attack. 
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Figure 3.3: Velocity and forces on a blade element. Wind Energy Handbook, Burton 
et al. Wiley 2001 [4]. 


3.3 Blade Element Momentum theory 
3.3.1 Drag and Lift forces 

Once the total velocity W of air particles past the single blade element is known, 
see equation (3.15), it becomes relevant to compute the aerodynamic forces acting 


on it. To this aim, the aerodynamic properties of the blades, namely the drag and 
lift coefficients Cp and Cz, respectively, lead us to write 


1 

ôD = PW eC dr (3.18) 
1 

ôL = 3PW°cCrdr (3.19) 


where ôr denotes an elemental ring belonging to the rotor plane and c is the blade 
chord. See figure 3.4. 


chord length 


Figure 3.4: Typical geometry of NACA airfoil. 


Together with the chord length c, the main geometric parameters which define 
an airfoil are shown in figure 3.4, where f is the maximum camber, x} is the position 
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of the maximum camber, d is the maximum airfoil thickness, xq is the position of 
maximum thickness, ry is the nose radius and finally yo (x) and y, (x) are the 
coordinates of the upper and lower side contours. 

As we are mainly interested in the forces normal and tangential to the rotor 
plane, the following projections are considered 


P,, = 6Lcos¢ + 6Dsing (3.20) 
P, = 6Lsin¢ — dD cos ọ (3.21) 


which can be normalized with respect to 1/20W?c to obtain 


Pn = dl cos d + ddsing (3.22) 
pi = dl sind — ôd cos (3.23) 
where 

Pn 

Pn = (3.24) 
4pW?c 

P; 
= — 3.25 
Pt i pW2c ( ) 


Projecting also the drag and lift coefficients we have 


Ch = Cr cos ọ + Cp sind (3.26) 
C, = Cr sing — Cp cos ġ (3.27) 


where Cn and C; are such as 


Il 


1 
Pn 59Cn Ww? (3.28) 


1 
pecs w? (3.29) 


II 


Pt 


Now, if N, denotes the number of blades, we obtain that the thrust and torque 
for an elemental volume of thickness ôr are, respectively 


hoo Wea 
ôF = 5 N n cOnòr (3.30) 
o1 W (1—a)Qr(1+a') 
ôM = 5PNe vasi cCyr or (3.31) 


By equating the two equations above, and by introducing the solidity o = 
crN,/ (27r), the system for a and a’ becomes 


1 
aCn + 1 
di (3.33) 
Asingcosd _ 1 


oC 
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Algorithm 1: Blade Element Momentum algorithm 
Data: Uso, 6, Q, r 
Result: ôF, 6M 
initialization (a = 0, a’ = 0); 
while toll (a,a’) <= toll* do 
Compute flow angle @; 
Read off Cz (a) and Cp (a) from tables; 
Compute Cn and C; from equations 3.26 and 3.27; 
Calculate a and a’ from equations 3.32 and 3.33; 
end 
Compute the local loads on the segment of the blades; 


The above system concludes the so called Blade Element Momentum (BEM) 
model which can be translated into a computer code by the following algorithm 1 
[34]. 

To sum up, the BEM theory just couples the momentum equations with the local 
element aerodynamics balance equations, involving drag and lift forces, to solve the 
complete system. 

The BEM method requires two corrections in order to account for both the real 
number of blades and the case when the momentum theory is no longer allowed. 
Finally, the corrected BEM theory can be satisfactory used to compute the loads 
on the rotor as well as the annual energy production. 

The aerodynamic solver used by the structural aero-elastic simulator which will 
be used in this thesis implements also the so called Generalized Dynamic Wake 
(GDW) model to compute the aerodynamic forces. This theory is also known as the 
Method of Acceleration Potential. It is an approach originally developed for heli- 
copter industry and exhibits some advantages with respect to the above discussed 
BEM. It allows a more general distribution of pressure across the rotor plane than 
the BEM. The GDW model is based an a solution of Laplace’equation for potential 
problems. Details about the GDW theory are available in [35], [4]. The aerodynamic 
solver used in the next simulations will use both the BEM and the GDW theories 
depending upon the current wind speed. 


3.4 Wind model 


Forces acting on offshore wind turbines mostly stem from the aerodynamics of 
the rotor, the offshore environment (waves, tides, currents, ice, ect.), the gravita- 
tional and inertial loads. Modern turbine rotor blades are getting larger and larger 
and this gives rise to an increase of the dimensions of all others structural compo- 
nents (e.g. the tower height). Larger dimensions cause, in turn, an augmentation of 
dead weight (gravitational loading), inertial forces and, last but not least, the effects 
of unsteady turbulent wind becomes more and more evident. A sketch of a three 
bladed onshore wind turbine in a full field turbulent wind is given in figure 3.5. 

In this thesis, as already pointed out in Chapter 1, an external condition—based 
extreme response analysis of offshore wind turbines is carried out. This means that 
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Figure 3.5: Effect of a turbulent wind speed distribution over the swept rotor area 
on the upwind velocity of the rotating rotor blades. Figure from [1]. 


extreme wind conditions are assumed according to recognized international stan- 
dards, such as IEC61400-1 and -3. 


Since the external conditions are dependent on the site, in IEC61400-1 wind 
turbines are divided in four classes. Each one is characterized by a reference wind 
speed U,e; and a turbulence intensity factor Ief. These wind classes have the intent 
to cover most of the onshore applications. While, a special wind class, referred to 
as “S” is devoted to offshore applications. No prescription are made on this special 
class where all parameters are specified by the designer. However, for offshore wind 
turbines, to design the Rotor—Nacelle Assembly, the definition of wind classes as 
in IEC61400-1 remains valid. Due to a lack of data, for the sake of simplicity, in the 


present study the parameters defining the class “S” are always chosen like those for 
wind class III. 
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Much attention must be paid to the reference wind speed, in fact if a turbine 
belongs to a specific wind turbine class with a reference wind speed Upes, it is 
designed to withstand climates for which the extreme 10 min average wind speed 
with a 50-year recurrence period at turbine hub height is lower than or equal to 
Uref- 

For the sake of simplicity, IEC61400-1 groups external condition in the so called 
Design Load Cases (DLCs) which are defined by combining: 


e normal design situations and appropriate normal or extreme external condi- 
tions; 


e fault design situations and appropriate external conditions; 


e transportation, installation and maintenance design situations and appropri- 
ate external conditions. 


It is worth pointing out that all structural and mechanical components are re- 
quired to resist both the ultimate and the fatigue loading conditions. While the 
design of tower and foundation is governed by the ultimate load cases, the design 
driver for the rotor blades is usually contemplated by the fatigue load cases. 

To each DLC it is assigned a specific type of analysis denoted by U (Ultimate), 
F (Fatigue). Ultimate analysis are additionally distinguished in Normal (N) or Ab- 
normal (A) and partial safety factors are then assigned accordingly?. 


3.4.1 Extreme turbulent wind speed model EWM 


As already pointed out, in this work the focus is on the extreme wind conditions, 
where the word “extreme” is referred to all those events with a 50-year return period. 
Offshore wind turbine must be designed to safely withstand wind conditions having 
intensity defined by such a return period. 

The randomness of the wind is taken into account by adopting the appropriate 
turbulence model. Among the two recommended by IEC61400-1 3rd edition, the 
Kaimal model is here adopted. The single-sided velocity spectra for the three wind 
components, k = u,v, w* is given as follows 


402 Ly /Unub 
(1 + 6fLx/Unun)5!* 


Sk (F) (3.34) 


where f is the frequency, ox is the standard deviation of the k-th velocity component 
and Lx is the integral scale parameter. See table 3.1 
where 


A= 0.7z z< 60m 
1) 42m 2>60m 


In this case the 10 min average wind speeds as functions of the elevation above 
the still water level, with 50-year and 1-year return period, respectively, are as- 


3As 1.35 for N, 1.1 for A situations. All fatigue design situations assume 1.0 as partial safety 
factor. 
4u is the longitudinal direction, v lateral and w vertical. 
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Velocity components 
| u v w 
Ou 0.80, 0.50, 
8.1A, 2.7, 0.66A; 


Ok 
Lk 


Table 3.1: Turbulence spectral parameter for the Kaimal model. 


sumed as follows 


r (3.35) 


Uy (2) = 0.8U re (3.36) 


The standard deviation in the longitudinal direction o, for the turbulent extreme 
wind speed model is 
Ou = 0.11U pub (3.37) 


where Unub is the wind velocity at hub height. 

Note that for specific design load conditions when the turbine is parked or stand- 
still (e.g. DLC 6.1a of IEC61400-3 [2]) the turbulent extreme wind model prescribed 
by the ICE61400-1 [36], see section 6.3.2.1, assumes that the turbulence standard 
deviation does not depend on the surface roughness. However, in general, this is not 
true and an appropriate model accounting for the effect of the sea surface roughness 
on the turbulence has to be adopted?. 


5 Apart from the extreme wind model, for other design load conditions, according to [2], the 
turbulence standard deviation, whenever there are no site data available, is related to the sea 
surface roughness as follows 


+ 1.28 x 1.45 x Is (3.38) 


where I15 is the average value of the hub height turbulence intensity determined at hub height 

wind speed Up, = 15 m/s. Offshore wind turbines are considered in the wind class “S” for which 

specific data regarding the wind and turbulence models have to be provided, see Annex A of [2]. 
While, the roughness length zo should be found by solving the following nonlinear equation 


2 
A U 
a _KUhub (3.39) 
g In ( Zhub ) 
Z 
where g is the acceleration due to gravity, k = 0.4 is the von Karman constant, Ac is the Charnock’s 


constant whose recommended values are 0.011 for open sea conditions and 0.034 for near-costal 
areas. See [2] for further details. 


Chapter4 
Hydrodynamic model 


This chapter describes the numerical model developed in this thesis for fully nonlinear 
water waves simulations. In addition to the aerodynamic loads discussed in the previous 
chapter, the offshore environment provides additional significant loading actions such as 
wave loads, current loads, ice impacts and tides. The numerical model here presented aims 
at the systematical inclusion (see next chapter) into the design procedures of the effects 
stemming from extreme waves breaking against offshore structures. 

The first part is devoted to a brief review of the standard waves descriptions in order 
to adequately prepare the background for the new impact wave model. 


4.1 Waves description 


Sea waves are traditionally described by both a deterministic approach and a 
probabilistic model with respective advantages and disadvantage. For instance, the 
spectral approach permits to describe a random sea, but has the drawback that 
only linear waves can be represented. As a consequence, only forces stemming from 
linear wave theories can be derived. 

On the contrary, some nonlinearity magnitude can be taken into account when a 
deterministic single-harmonic wave is used. In other words, the two commonly used 
approaches implement respectively either 


e Regular nonlinear waves, or 
e Irregular linear waves; 


Fortunately, in most cases the nature of oceans can be very well described by 
the superposition of linear regular waves, and this makes the spectral approach an 
extraordinarily effective tool, in fact especially for long term loading condition, this 
approach fits fairly well the nature of the actions. 

Contrarily, for ultimate failure conditions more representative models are re- 
quired in order to capture the fully nonlinear contribution due to extreme (possibly 
breaking) waves. 

Nowadays capabilities of modern computers permit to simulate fully nonlinear 
waves without penalizing the total simulation time, thus, whenever the fully non- 
linear behaviour of waves plays a dominant role in designing offshore structures, 
it seems to be opportune to adopt a direct numerical solution of the governing 
equation without introducing any a priori hypotheses!. 


1Note that “a priori hypotheses” refers to the magnitude of nonlinearity, rather than to the 
assumptions on the physical characteristic of the fluid and flow type. 
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Analytical theories, such as Stokes-2nd and 5th order, fall in the weakly non- 
linear group. In fact, what really distinguishes weak nonlinearity to the high (or 
strong) nonlinearity is essentially the asymmetry with respect to the vertical axis 
[31]. In the weakly nonlinear case we only deal with asymmetry with respect to the 
horizontal axis: crests become narrower and the troughs wider. 

Figure 4.1(a) shows the two traditional approaches used to describe ocean waves: 
deterministic and stochastic. The advantage of the deterministic description is the 
capability of simulating nonlinear waves both through analytical and numerical 
tools. The main drawback however is that the real nature of oceans is not de- 
terministic. On the other hand, as already pointed out, the stochastic approach 
permits to capture the real nature of the sea but only superimposing linear regular 
waves. To overcome the disadvantages of both methodologies, first a deterministic 
fully nonlinear Boundary Element Method-based simulator is developed, then it is 
linked with the wave energy spectra in order to account for the stochastic nature of 
the reality, see figure 4.1(b). 

In this thesis the link between the probabilistic and deterministic approaches is 
used to get the advantages from both by virtue of their complementarity. In fact, 
in a random wave series there are some larger components, characterized by special 
energy contents, which have to be be captured and described separately because the 
importance of their strong nonlinearity is dominant with respect to the whole sea 
surface representation. These are the cases when exceptional events, such as rouge 
waves, occur. An example is shown in figure 4.2. 

Figure 4.3, The Great Wave, shows how the danger represented by such extreme 
waves was also well impressed into artists’ imaginary. 

In the next section the fully nonlinear numerical wave model is described. Later 
on, in the next chapter, the numerical model is embedded into the stochastic envi- 
ronment. 


4.1.1 Deterministic representation 


There are several wave theories that provide a deterministic description of water 
waves. For most of them, three parameters are fundamental, while all the others 
can be derived from these. They are the period T, the wave height H and the water 
depth d. 

The most used wave theories are shortly listed below: 


Linear (or Airy) wave theory: the most important, useful and applied theory. 
It is also known as small amplitude wave theory and it is based on a strong 
linearization which makes the theory suitable for the probabilistic spectral 
representation of random seas. The analytical solution for the velocity poten- 
tial and all the kinematic quantities is found by dropping all the second order 
terms in the dynamic and kinematic boundary conditions at the free surface. 


It is always useful, especially for the scopes of Chapter 5, to have at hand 
all the kinematic and dynamic equations governing the propagation of a reg- 
ular sinusoidal wave. For this reason, the essential formulas are listed in ap- 
pendix A. 


2nd and 5th order Stokes theories: these theories are also known as finite am- 
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(b) Proposed approach. 


Figure 4.1: Traditional scheme and proposed analysis approach adopted for describ- 
ing ocean waves. 


plitude wave theories. In fact, they employ a perturbation parameter called 
steepness € = ka, where k is the wave number and a the wave amplitude, 
which permits to describe steeper waves. 


Cnoidal theory: the above Stokes theories have some restrictions on the applica- 
bility in shallow waters. The cnoidal theory supplies a proper description for 
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(a) Deep water rogue wave. (b) Shallow water freak wave. 


Figure 4.2: Examples of freak waves. 


Figure 4.3: The Great Wave by Katsushika Hokusai, 1760-1849. 


finite amplitude long waves in shallow waters. A cnoidal wave has a typical 
shape consisting of sharper crest separated by wide troughs. 


The importance of understanding the difference between several wave theories 
comes up if we look at the applicability diagram in figure 4.4. 

A complete description of all the theories mentioned above is available in many 
books [37], [5], among others. 


4.1.2 Probabilistic representation 


As already mentioned, to design offshore structures both the stochastic and the 
deterministic approaches are necessary. To describe the random nature of ocean 
waves the only tool is represented by the wave energy spectrum. Energy density 
spectra represent the energy content of an ocean wave and how it varies on a certain 
range of frequencies. 
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Figure 4.4: Wave theories applicability, from [5]. 


It is not our aim to repeat what is available in a very wide range of textbooks, but 
to have at hand some important concepts we just mention the two most important 
wave spectra which have been implemented in the numerical model discussed in the 
next chapter. 


Standard wave spectra 


After carrying out an environmental analysis, the key parameters which describe 
the sea severity can be established. They are: the wind velocity, the significant wave 
height, the mean zero—crossing wave period. Then, depending upon the type of sea 
to be simulated, the two most used wave spectra are [32], [11], [31], [38], [2], [7]: 


Pierson—Moskowitz: The Pierson-Moskowitz spectrum was developed departing 
from data obtained in the North Atlantic in condition of fully developed sea. 
This is a single-parameter formulation, indeed it depends only on the wind 
velocity U measured at 19.5m above the sea water level. 


8.10 g? (#4) 
.10 g a 
Sin (w) = 103 we oP 0.032 J (4.1) 
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where H, is the significative wave height related to the wind velocity U as 
H, = 0.21%. 


JONSWAP: The Joint North Sea Wave Project spectrum is based on an extensive 
measurement campaign carried out between the 1968 and 1969 along a line of 
about 160 km in the area of Sylt Island, see [11] and its references. Contrarily 
to the Pierson-Moskowitz formulation, this spectrum is suitable for wind— 
generated fetch-limited seas. Therefore it requires as input both the wind 
velocity and the fetch length. The formulation is given as follows 


2 cal] 


_ nl oxy | (#0) exp|- Seat 
Sm (w) = a |- (2) | (4.2) 


where 


e w, is the peak circular frequency. According to [11] it is related to the 
fetch F and wind velocity U by the following relation 


Wp = Tiap” (4.3) 


e 7 is the peak-shape parameter, it represents the ratio between the max- 
imum spectral energy density and the maximum of the corresponding 
Pierson—Moskowitz spectrum 


e a is a parameter related to the sea generation conditions, indeed it de- 
pends on the fetch F as follows 
a = 0.076170”? (4.4) 


e c = 0.07 for w < wp and o = 0.09 for o > wp 
o rp = A is the nondimensional fetch length 
e U is the mean wind speed 

e F is the fetch length 


e gis the gravity acceleration 


It is always useful, however, to obtain the JONSWAP spectrum in terms of the sea 
severity H,. In agreement with [11] and [12], we will assume the following relation 


Se ieee eae (4.5) 


where k depends on the peak shape parameter y, as shown in table 4.1, F is the 
fetch in km and H, must be expressed in m. 

The line in red in table 4.1 highlights the values for k associated with the mean 
value of y. 

The value y varies approximately from 1 to 5 randomly, it is usually normal 
distributed with mean value of 3.30. However, according to [2] and [10], the following 
formulation is adopted: 


5 o for vir < < 
yY=$ exp (5.75 <1, He) for 3.6 < < Jir < 
1 for - > 5 
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k y 
96.2 1.75 
88.3 2.64 
83.7 3.30 
80.1 3.96 
76.4 4.85 


Table 4.1: k — y relation, [11] and [12]. 


It is convenient sometime to get the JONSWAP spectrum as function of the two 
independent variables Hs and Tp. In this way we do not use neither table 4.1 nor 
the fetch F which is an information already contained in 7, and H,. Indeed, given 
these two inputs, a is readily obtained as follows [39] 


1 2 
1 \1. [A 1 H 
A ESE ay | ay ue ee 4. 
p (sas) Agi en l 9 ) (he) 


In conclusion, the formulation given in equation (4.2) would require as indepen- 
dent variables U and F. In other design situations Hs, and T, are assumed the two 
independent design parameters. From them, of course, it is possible to get back first 
a via (4.6), then F through (4.4) and finally U. 


4.2 Fully nonlinear potential flow water waves 


Most of the current approaches in designing offshore structures successfully 
adopt the linear wave theory, nevertheless in some design conditions nonlinear effects 
cannot be neglected, especially when the goal is evaluating the structural safety. 

Although modern computer simulations of offshore wind turbines have made 
formidable progresses, the integration of a fully nonlinear numerical solution of 
gravity waves into the more general multi-physics framework characterizing the 
design of offshore wind turbines seems to be not yet a common practice. 

One of the first contribution addressing the numerical simulation of nonlinear 
water waves was due to Longuet-Higgins and Cokelet in [40] who introduced for 
the first time the Mixed Eulerian-Lagrangian approach to describe such a free 
surface problem. Subsequently, in [41] it was proposed a new and time-effective 
procedure to integrate in time the dynamic and kinematic boundary conditions 
on the free surface. The advantage of this time-stepping procedure, indeed, lies in 
solving different Laplace’s equation at each fixed time step making use of the same 
system matrices. 

The numerical solution of the Boundary Value Problem (BVP) and the conse- 
quent time integration do not introduce any approximation neither on the velocity 
potential nor on the dimension of typical wave parameters such as wave height, 
wave period, water depth and wave length. Excellent results of this approach in the 
description of fully nonlinear water waves have been achieved in [42], [43], [44] and 
[45], just to mention a few. 

Nowadays the implementation of this computational tool deserves attention for 
several reasons: firstly we detected a lack of fully nonlinear models in investigat- 
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ing structure-fluid interaction for offshore wind turbines substructures. Secondly 
because, especially in those seas where the water depth is fairly high, floating struc- 
tures are desirable and in this case methods hereafter illustrated are quite promising 
[46], [47] and [48]. Moreover, as for long term loading actions the sea state needs 
to be described by a superposition of nonlinear waves, interesting developments are 
also in the direction proposed in [49]. 


4.2.1 On the validity of the potential flow model to describe breaking waves 


The potential flow theory is regularly used to model water waves. The existence 
of a velocity potential stems from the hypothesis that the fluid is irrotational. In 
the model used in this work the condition of non-rotational flow is guaranteed 
throughout the evolution of the wave up to the time instant when the water jet 
re-enters into the free surface. At the re-entry time instant, in fact, the domain 
becomes multi-connected and Kelvin’s theorem does not hold anymore. 


Kelvin’s theorem states that for an ideal fluid subjected to conservative body 
forces, the circulation about any closed material contour moving with the fluid is 
constant in time, see e.g. [50]. 


So that, given a zero initial circulation, we are sure that it remains zero up to the 
re-entry. Moreover, Stokes’ theorem assures that the condition of zero circulation 
on the closed contour is equivalent to the condition of irrotaional flow. 


In addition to the theoretical justifications, experimental evidences (see e.g. 
[6, 51]) confirm that the irrotational model is valid to capture the evolution of 
plunging breakers (up to the re-entry). 


The above argumentations lead to the conclusion that no relevant discrepancies 
from the real flow situation are caused by the assumption of a potential model. 


As already mentioned, in addition to the hypotheses of incompressibility of the 
fluid and irrotational flow discussed above, the present model is based on the as- 
sumption of non-viscous fluid. It has been demonstrated that the effects of the vis- 
cosity on the kinematics of water waves is negligible for the cases we are interested 
in. In general, the viscosity might influence by acting through three mechanisms: 
1) viscous effects at the free surface; 2) viscous effects inside the fluid; 3) viscous 
effects at the sea bottom. Concerning mechanisms 1) and 2), viscosity causes the 
progressive attenuation of the gravity waves, but this happens for time periods much 
longer than those typical for the propagation of gravity waves. In case 3), the effect 
could be taken in some consideration only in the case of very small water depth. 
Our simulations are all in intermediate water conditions, thus the effects of viscosity 
can be neglected without introducing significant errors. Further details on the effect 
of viscosity can be found in [52]. 


On the contrary, the viscosity plays a different role in the loading model. In fact, 
a non-viscous fluid would lead to no drag contribution in the hydrodynamic force 
exerted on the monopile. To avoid this, Morison’s equation provides the drag term 
which accounts for the viscosity of the fluid. 
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4.2.2 Governing equations 


From Euler’s equations valid for an incompressible and inviscid fluid, the addi- 
tional hypothesis of irrotational flow allows the description of fully nonlinear water 
waves by means of a potential model. At a fixed time t, the velocity field for each 
point belonging to the spatial domain Q (t) can be expressed by means of a velocity 
potential ¢ (p, t) as follows 


u (p,t)=Vọ(p,t) Vo Ee Q(t) (4.7) 


Equation (4.7), together with the mass conservation (div (v) = 0) leads to the 
following Laplace’s equation 


V°o(p,t) =0 Vp Ee Q(t) (4.8) 


The domain Q (t) is assumed bounded by four boundaries: T; (t), Ts (t), Do (t), 
I’; (t) being the inflow wall, the bottom, the outflow wall and the free surface, 
respectively. See figure 4.5. 


I) 


T(t) 


I 


Figure 4.5: Two-dimensional domain of the potential problem. 


Let n be the unit outward normal vector, the normal component of the velocity 
field (the flux) stems form (4.7), as follows 


Equation (4.9) is used to provide Neumann boundary conditions on the bound- 
ary T; (t), Ty (t), To (t). The value to be assigned to the normal derivative of the 
velocity potential depends on the type of water wave to be simulated. 

An inertial coordinate system is fixed with the x-axis along the Still Water Level 
(SWL) and the y-axis vertical and upwardly oriented. See figure 4.5. At a fixed time, 
points p € I ș are tracked by a Lagrangian position vector 7 (p) = p—o = x;èx+yfèy. 
Where e, and €, denote the unit normal basis. 
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Boundary conditions 


On the free surface T ș (t) the dynamic boundary condition, which represents the 
Bernoulli equation for unsteady flow, reads as follows 


Do (pt) _ _ Pa 
Dt Pw 


gys + 3700) Volt) Versi) (410) 


where yy is the free surface instantaneous elevation, often also denoted by n, pa is 
the atmospheric pressure, which can be neglected, and p,, the water density. 

Moreover, the profile of the free surface is governed by the following kinematic 
boundary condition which stems directly from equation (4.7) 


DIP) 2. 
PFO) _ 5 (p,)=Volp.t) Yoel y(t (4.11) 
The above free surface kinematic boundary condition can be written in compo- 


nents as follows 


Dx; i 
DI = Uf (4.12) 
Dyş y 


To solve the time-depending potential problem (4.8), initial conditions must be 
also assigned. Namely, at the beginning of the simulation the potential along the 
free surface and its geometry have to be known. 


4.2.3 Time integration scheme 


There are several methods to integrate the dynamic and kinematic boundary 
conditions on the free surface. A review of the possible approaches is proposed in 
[53], [54]. One of the most effective time integration scheme for this type of problems 
was first proposed by Dold and Peregrine [41]. It is based on Taylor series truncated 
at a certain order which permits to approximate, and consequently to update, both 
the free surface profile and the velocity potential from the current time t to the 
subsequent time step t + dt. Namely, the series are 


Dr (p,t) , , 1 D°r(p,t) 


F(p,t+dt)=F(p,t)+ pet po dt? + o (dt?) (4.14) 
2 
ġ (p, t + dt) = @ (p, t) 4 PotD a 4 i OD ap + o (at?) (4.15) 


Figure 4.6 shows the Lagrangian updating of the free surface water particles 
position. 

The above series are truncated at the second order. It seems that such a choice 
is optimal considering the numerical effort necessary to compute the third order 
coefficients [55]. The procedure to get up to the second order coefficients of the 
above series is rather simple and has some remarkable advantages with respect 
to other time-stepping schemes. The following two paragraphes describe how to 
proceed. 
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Figure 4.6: Lagrangian updating of the free surface particles position. 


First order Lagrangian derivatives 


No special effort is required to compute the first-order coefficients in equations 
(4.14) and (4.15). Indeed, the first-order total derivative of the position vector 7, 
i.e. the velocity components vy and v$, can be computed by accounting for both 
the contributions of normal and tangential components 


vp = vpn” + vpt” (4.16) 
vy = vn” + vpt” (4.17) 


where, in turn, vp stems directly form the solution of the BVP, while the tangential 
component vh is obtained by differentiating the shape functions. See section B.2 of 
appendix B. 

In equations (4.16) and (4.17), n”, n” and t”, t¥ denote the Cartesian components 
of the normal and tangential unit vectors, respectively. 

In addition to this, the first-order coefficient for the integration of the velocity 
potential, that is its total derivative, see equation (4.15), can be directly computed 
by means of the dynamic boundary condition (4.10) 


Dé (p, t) "Ea I 1 x2 | y2 
Di MF +5 (v H vy) (4.18) 


Second order Lagrangian derivatives 


A bit more articulated procedure is invoked to compute the second-order coef- 
ficients involved in equations (4.14) and (4.15), [45]. Let us start from the particle 
acceleration 

Dr Dü _ dv 


DE = DI pp + VD ©) (4.19) 


By means of equation (4.7), the temporal derivative a can be written as follows 
dv ; 
— =V 4.20 
Två (4.20) 
where it has been set ġ — de Therefore, the particle acceleration becomes 
Dt ; 
2 = Vo + (V5) (7) (4.21) 
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To solve the above problem it is necessary to know ò. By differentiating equation 
(4.8) with respect to time, it is immediate to realize that for a fixed time, @ satisfies 
another potential problem stated as follows 


V26(p)=0 Vpea (4.22) 
oe Do = (v7? + wp) Vpely (4.23) 
ùp = Si Yp ETs (4.24) 
i? = Vpell (4.25) 
oP = sae Vper, (4.26) 


Note that the right-hand side of equation (4.23) is completely known since the 
total derivative of the potential stems directly form the dynamic boundary condition 
and the velocity field on T has already been computed by using equations (4.16) 
and (4.17). By solving this second BVP, with the same distribution of Dirichlet and 
Neumann boundary conditions, we get the unknowns OF. Furthermore, by numerical 
differentiation it is possible to calculate also oF and thus it is straightforward to 


compute the gradient of O) as follows 


dd » nN Xx ot 4X 

da = Vf = Vf ST ùrt (4.27) 
OP ci o 

y =ù} = Un + bpt” (4.28) 


The above approach is particularly effective because the system matrix needs to 
be computed only once at each time-step as the geometry of the domain £ is the 
same for both BVPs. Namely, this permits to use again the same matrices H and 
G, see equation (B.28) in appendix B. The difference between the two BVPs lies 
only in the boundary condition value (not type!), and so no new integration over 
the boundary elements is required. 

The second-order coefficient for Taylor series (4.14) still need the last contri- 
bution (Vv) (v). The procedure to compute this involves nothing but some basic 
vector calculus where the irrotational hypothesis of the flow and the continuity equa- 
tion are used [45]. However, details on the calculation of this term are discussed in 
section B.3 of appendix B. 

Concerning the second-order coefficient a necessary to update the velocity 
potential at time t+dt, it is required to differentiate the dynamic boundary condition 
with respect to time as follows 


D D 1o y a DoF „DoF 
De T 5 ( oi 3%) PA DE PSE ee) 


4.2.4 Method of solution 


As already mentioned, the fully nonlinear potential flow initial-boundary-value 
problem is solved by using the Mixed Eulerian Lagrangian (MEL) approach which 
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consists in a repeated two-step procedure. See [56] for a wide literature survey. 

After the notions introduced in the previous sections, the MEL approach can 
now be better highlighted. 

For the sake of brevity let Ty =T; UT, UT, denote the region of the boundary 
where Neumann boundary conditions are assigned. Note that the dependency of 
the boundary curve I on time has been suppressed. This abuse in notation is here 
justified by the fact that at each time step the problem is regarded as a steady case. 

The two steps are: 


1. Eulerian step: at a fixed time t, the free surfac I°; and the velocity poten- 
tial ¢(p € T;) on it are known. The flux Vọ (pe T;)- on Ty and Ty it- 
self are also known. Two steady Laplace’s equations are solved for the fluxes 
Vé(per;)- f and Vò (p € T): ^ñ on the free surface and for the velocity 
potentials @ and $ on Ty. 


2. Lagrangian step: making use of equations (4.14) and (4.15), the velocity po- 
tential and the free profile are updated in time providing boundary conditions 
for the next Eulerian step. See figure 4.6. 


At the second step, when all the unknowns have been found, each particle of the 
boundary of the Eulerian frame is updated in a Lagrangian manner. 

It is worth pointing out that the particular type of time-integration adopted 
requires that at each time step the number of steady Laplace’s equations to be 
solved depends on the order of Taylor’s series, so that, in accordance with equations 
(4.14) and (4.15), in the current code two BVPs are solved at each time-step. 

The steady solution at each time-step is achieved by discretizing the BVP by 
using the direct Boundary Elements Method. See appendix B.1. 

Note that in the above procedure the analytic linear solution (Airy wave theory 
for irregular waves) plays a twofold role: it is used to initialize the solver by providing 
the free surface elevation and the velocity potential on the free surface, and it is 
used to provide Neumann boundary conditions (involved in the Eulerian step) on 
the upstream and downstream walls during the simulation. Moreover, the transition 
from the linear solution to the fully nonlinear one is made by using a ramp function 
which is required to be long not more than ten times one boundary element length. 
In Chapter 5, in particular in section 5.4.3, such a domain decomposition strategy 
will be extensively discussed. 


4.2.5 Smoothing and regridding 


According to [57| in our model the two most important and typical numeri- 
cal instabilities occurred: (i) strong instability due to a too large time-step; (ii) 
steep wave instability, also known as “sawtooth” instability. The Boundary Element 
Method used to discretize Laplace’s equation implements second order elements 
along all the four boundaries. Hence, the sawtooth instabilities look like a bit dif- 
ferent to the classical one shown in literature. Figure 4.7 gives an example of what 
this instability causes. 

While the strong instability can be fixed by setting a proper time-step, to remove 
the sawtooth behavior of the free surface a smoothing procedure is necessary. The 
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Figure 4.7: Typical sawtooth instability affecting the free surface in the case of a 
steep wave generated by a piston wavemaker. 


are many possibilities to smooth the free surface, here the one proposed in [40] has 
been chosen. Every a variable number of time steps the free surface and the velocity 
potential are smoothed by the following 5—point smoothing formula 


1 
J 16 


smt ana 


(-fi-2+4f;-1+10f;+4fj41— f;+2) (4.30) 


where ie is the smoothed value in the j-th node of the function f which can be 
$, yf and zp. 


The above smoothing formula is valid for equally spaced points and does not 
apply to the first (and last) two nodes of the free surface. When nodes tend to 
gather in the large velocity gradient regions, the equally spaced condition is no 
longer sufficiently accurate, hence the smoothing subroutine also implements the 
generalized form of the 5-point scheme in accordance with [58]. 


When dealing with overturning waves, however, the formation of a water jet 
causes nodes to concentrate near the cusp where very high velocity fields are asso- 
ciated with high curvatures. In this situation a more refined mesh could be needed 
but it must also be avoided that nodes undergo displacement so large to step over 
the neighbor particles. A regridding subroutine has been thus implemented with 
the aim of avoiding the latter inconvenient. When necessary, the regridding is also 
used to refine the free surface discretization as it allows the augmentation of the 
boundary elements on selected boundary patches. The regridding makes use of cubic 
splines. 


Smoothing and regridding have to be used together carefully. For both there are 
advantages and drawbacks which can lead to inaccurate solutions. 
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4.3 Validation 


4.3.1 Periodic waves 


The first simplest case used to test the numerical model regards a period wave. 
It will result clear later on the importance of demanding excellent performance to 
the code in simulating linear periodic waves. In fact, when the goal will be the 
simulation of an irregular sea, it will be just matter of the number of harmonics 
of which the velocity potential and the other kinematic entities are made up, but 
the global framework of the numerical solver is absolutely identical to this case. 
However, in this specific case, to check also the capability of the code of capturing 
nonlinear effect, a periodic 2nd-order Stokes wave has been simulated. Care needs 
to be devoted to the inflow and outflow boundaries where periodicity guarantees 
no reflection. Although the analytical solution for such a problem is spread into an 
unbounded number of textbooks, e.g. [37], here the initial and boundary conditions 
are recalled and adapted to the present numerical scheme. 


Boundary conditions on Ty 


The boundary condition on the bottom is the no flux condition, so that 


ve =Vb-n=0 Vpery (4.31) 


Boundary conditions on T; and To 


The velocity along the inflow and outflow boundary are consistent with Stokes 
2nd-order theory. They are 


i uti cos (k (ai — et) + 
? cosh (kd) 
A =) o kd) (2 + cos (2kd)) cos (2k (x; — ct)) (4.32) 
a Hs en s (k (£o — ct)) + 


h 
2(4 =) ne aH | 2 + cos (2kd)) cos (2k (£o —ct)) (4.33) 


where s = d + y. 
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Initial conditions 


The initial conditions are assigned to the free surface as follows 


Hg cosh (ks) 
~ 2w cosh (kd) 
3 H?n? cosh (2ks) 
8 KTL sinh? (kd) 


@ (p, 0) sin (kx) + 


sin (2x) Vp €T (0) (4.34) 


A 
yf (2,0) = y COS (ka) + 


nH? cosh (kd) 
8L sinb? (kd) 


(2 + cos (2kd)) cos (2kx) Vp € Ty (0) (4.35) 


where yy; (often denoted also by 7) is the elevation of the free surface with respect 
to the still water level. 


4.3.2 Solitary wave 


The second numerical experiments to test our code addresses the propagation 
of a solitary wave as described in [59]. Actually, several kinds of solitary waves have 
been simulated, e.g. [55] among others, and for all cases excellent results have been 
attained. Here, for the sake of brevity, only the solitary wave simulated in [59] is 
shown. 


Boundary conditions on T;, To and Ty 


To generate a solitary wave we give some specific initial conditions on the free 
surface, while all the remaining boundaries are assumed impervious and kept at rest 
by setting the following Neumann boundary conditions 


vp (Vi, t) = v; (Ti, t) = vp (Tot) =0 Vie [t:i tr] (4.36) 
where t; and t; are the initial and final instants of the simulation. 
Initial conditions 


The motion is generated by initial conditions obtained by the following second- 
order analytical solution setting t = 0 


$t) =s 
tanh (u — kàt) V4/3H (1 — 5/4H)+ 

+ 1/Asech? (u — KAt)? tanh (u — kàt) \/3/4H3 (1 = 5/4H) (4.37) 

n (x,t) = H [1 — 3/4H tanh? (4 — KAt)] sech? (u — kAt) (4.38) 


[1 — 1/4H +1/3H tanh? (u — kat)’ x 
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where 
k = \/3/4H (1 — 5/4H) (4.39) 
A=1+1/2H — 3/20H° (4.40) 
u= ka (4.41) 


4.3.3 Piston wavemaker 


Here a fully nonlinear numerical wave tank is simulated. See [60], [61], [62], [63], 
[43], [42], [51], [46], [6], [64], [45], [44], [65], among others. 

Two fundamental cases have been tested: (i) when the paddle is moving gen- 
erating single harmonic wave and (ii) the motion of the piston is given by the 
superposition of several harmonics. 

The numerical piston starts moving at the initial time t; = 0 when the fluid lies 
in a state of rest with the free surface being horizontal. 


Boundary conditions on T, and VT, 


The outflow wall T, and the bottom T, are not time-depending and for both of 
them the no-flux condition is assigned as follows 

uy =Vo-n=0 VpeTy (4.42) 

ve =Vo-n=0 Vpero (4.43) 


Boundary conditions on the wavemaker 1; 


According to [51], [6], [64], in the case of single harmonic, the piston moves with 
the following general law 


£p (t) = i cos (wt) (4.44) 


so that the position of the piston at the initial instant is given by the negative 
semistroke —A/2. From equation (4.44), the upstream flux assumes the following 
expression 


Up (t) = du sin (wt) (4.45) 


Note that since the unit normal vector is always outwardly oriented from the 
domain Q (t), the Neumann boundary condition assumes the negative sign: v? (t) = 


—Up (t). 
Initial conditions 
The fluid in the water tank is at rest, so that the initial conditions on the free 


surface are assumed as follows 
¢(p,0) =0 Vp els (0) (4.46) 
f(x,0)=0 VpeET; (0) (4.47) 


where f is the elevation of the free surface with respect to the still water level, 
indeed we have yy (x+,t) = f (x,t). 
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4.3.4 Absorbing beach 


To reproduce a real wave tank, especially when simulating waves generated by 
piston wavemakers, it is necessary to avoid the reflection of the energy introduced 
into the system by the piston. This can be done by means of a damping layer which 
can be tuned in order do dissipate a certain amount of energy. See figure 4.8. 

Thus, the waves are gradually attenuated along the beach which has a fixed 
length. Along this sponge layer both the kinematic and dynamic free surface bound- 
ary conditions have been modified by introducing the dissipative term as suggested 
in [64] as follows 


Pel: 2 7 gys + 5VO(0, t)-Vé(p,t)-v(p)(d- de) Wp eT (t) (4.48) 


and for the free surface kinematic boundary condition (4.11), we have 


DF (p,t) 


pi Ut) = Velp t)-v(p)(T-7.) Vee s(t) VpeT y(t) (4.49) 


where ġe and Fe denote the velocity potential and the free surface particle position 
at the reference configuration, that is when no waves are being generated. 

The absorbing piston located downstream at I’, as described in [61] is not nec- 
essary for our purposes. 

Setting the function v (p) is crucial as its strength makes the filter more or less 
effective. The scope of absorbing the incident wave energy before it reaches the 
wall and thus being reflected depends on the parameters which tune the following 
expression 


v (£) = aw È (a — DI “ype [zo, x1] (4.50) 


where xo is the starting point of the beach and zı = zo + 27 is the channel length. 
To fully absorb a wave characterized by w and k as angular frequency and wave 
number respectively, the parameters a and 8 should both equal one. 
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Figure 4.8: Sketch of the numerical wave tank equipped with a sponge layer. 


4.4 Results 
4.4.1 Stokes 2nd-order 


The simulated wave is characterized by H = 0.2m, d = 2.5m, L = 6.28m, 
T = 2.019s. Where H denotes the wave height, d the water depth, L the wave 
length and T the period. The harmonic wave is simulated considering a numerical 
tank 4L = 25.12m long?. 

Figure 4.9 shows the comparison between analytical and numerical free surface 
evolution obtained with NE; = 60, NE, = 20, NE; = NE, = 10, dt = 0.05s. 
Where NE denotes the number of boundary elements on the a-boundary with 
a = i,b,o and i: inflow (upstream); b: bottom; o: outflow (downstream). Moreover, 
the smoothing scheme introduced above has been used. 


4.4.2 Solitary wave 


The parameters defining the solitary wave and the geometry of the numerical 
tank are 


e wave height H = 0.1m; 

e water depth d= 1 m; 

The number of elements on the free surface is NE; = 48, while on the upstream, 
downstream and on the bottom there are NE; = NE, = 8 and NE, = 24 elements, 


respectively. The time-step adopted to integrate the boundary conditions is dt = 
0.05s. 


Note that H/ (gT?) = 0.005 and d/ (gT?) = 0.0625. These values make the simulated wave 
falling in the Stokes second order theory applicability, [37]. 
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Figure 4.9: Analytical and numerical free wave propagation of a second order Stokes 
periodic wave. Free surfaces taken at t = 2s. 


To generate a solitary wave we assigned initial conditions from the analytical 
solution for t = 0. As said earlier, on the upstream, bottom and downstream walls 
the boundary conditions are defined by the impervious (no flux) condition. 

Figure 4.10 shows the excellent agreement between the numerical results and 
the analytical solution. In addition to this, a proof of the efficiency of the code is 
given by computing both the total mass of the system and the flux trough the free 
surface. In the former the mass conservation is always guaranteed and in the latter 
zero flux balance is also preserved (with 10E-7 approximation). 

The comparison presented in figure 4.10 has been intentionally limited to t = 
11, because just further this dimensionless instant the run up starts, as shown in 
figures 4.11(c) and 4.11(d), and the comparison becomes meaningless. 

A complete propagation is presented in figure 4.11 where also the run up is well 
visible. The simulation has been stopped at t = 30. 

Simulation in figure 4.11 has been run with the following parameters: NE}; = 64, 
NE, = 20; NE; = NE, = 10; dt = 0.05; d = 1. All the parameters have been made 
dimensionless by using water depth d for lengths and gd for time. 

As already mentioned, along with fundamental comparisons with analytical re- 
sults, some global quantities are also controlled to be in agreement with the theo- 
retical expected values, naturally within some tolerance. 

First, the total mass involved in the system has to be unchanged and this is 
shown in figure 4.12. Another useful check of the numerical reliability concerns 
the total flux balance. By using the divergence theorem, indeed, given Laplace’s 
equation on a domain with boundary T, then it is possible to prove that 


NE; 


96 ap ; G) 
: T = 2 f. X yr (s) dla = 0 (4.51) 


i k=1 


where x (s) is the k-th shape function and NE; the number of quadratic elements 
used to discretize the free surface. See appendix B for further details. Figure 4.13 
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Figure 4.10: Analytical and numerical free wave propagation of a solitary wave. Free 
surfaces taken at each 40 time steps (At = 40dt = 2, from t = 0 to t = 11. 


shows total flux balance during the propagation of the solitary wave. 
4.4.3 Piston wavemaker: regular wave 


The parameters concerning this simulation are itemized in the following. Note 
that all depends on the wave length L and the piston stroke S. The remaining 
parameters have been derived accordingly. The wave height in particular has been 
computed by invoking the linear transfer function. 


e tank length L; = 32 m; 

e water depth d = 1m; 

e wave length L = 8m; 

e piston stroke S = 2.0 x 1071 m 

e wave number k = 7.854 x 107}; 
e angular frequency w = 2.248 rad/s; 
e wave period T = 2.795 s 

e celerity c = 2.862 m/s; 

e wave height H = 1.559 x 107! m; 
e d/(gT?) = 1.305 x 1072; 

e H/(gT°)= 2.034 x 1073 


e initial simulation time t; = 0; 
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(a) Free surface at t = 0. 


BEM Free Surface at t = 15.0000 


(b) Free surface at t = 5. 


BEM Free Surface at t = 17.5000 


(c) Free surface at t = 10. 
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(g) Free surface at t = 22.5. 
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(i) Free surface at t = 30. 


Figure 4.11: Propagation and run-up on an vertical wall of a solitary wave. 
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Figure 4.12: Mass conservation during the propagation of a solitary wave. 
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Figure 4.13: Total flux through the boundary during the propagation and runup of 
a solitary wave. 


e final simulation time t; = T; 


Figures 4.14 to 4.18 show the evolution of a nonlinear regular wave generated by 
a piston wavemaker. They show the transient and the beginning of the stationary 
evolution. Even though linear theory is adopted for the transfer function, expected 
results in terms of wave length and wave height are excellently met. 

This simulation implements also the absorbing beach which is fundamental in 
these type of simulations in order to avoid reflection. When a regular wave is being 
simulated, setting the numerical beach parameters is straightforward as the sponge 
layer is calibrated to absorb the whole energy associated with that single harmonic. 
In this case we have in fact 


ea=1; 
e B=1; 
e % = 24m; 
e zı = 32m; 


where it results that the beach length just equals the wave length. 

The number of elements on the free surface is NE; = 60, while on the upstream 
(the paddle) and downstream boundaries we have NE; = NE, = 2. On the bottom 
there are NE, = 15 elements. The time-step adopted to integrate the boundary 
conditions is dt = 0.05s. 


4.4.4 Piston wavemaker: breaking wave 


In this paragraph we show the capability of the code of simulating breaking 
waves induced by wave-wave interaction. Contrarily to the previous cases, where 
the augmentation of steepness is well controlled by the stroke of the piston, in this 
case we are to face a real fully nonlinear phenomenon which is characterized by very 
rapid increase of velocity together with high curvature regions of the free surface. 
This scheme is numerically very sensitive and a high resolution is required both in 
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Figure 4.14: Propagation of a regular wave generated by a piston wavemaker. 


space and in time whenever the wave becomes unstable. The refinement in space 
is gained by implementing a local mesh refinement subroutine while the latter is 
achieved by halving the integration time step when each node moves more than a 
fixed percentage of the adjacent elements length. 

A deep water plunging breaking wave is here simulated according to [6], where 
an irregular wave field is generated by the following motion of the paddle 


72 
Up (t) = X. Un cos (wnt — On) (4.52) 


where vp is the velocity of the piston. The velocities U,, the angular frequencies wn 
and On the phases induce a breaking wave due to wave-wave interaction at 11.5m 
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Figure 4.15: Propagation of a regular wave generated by a piston wavemaker. 


approximately apart from the paddle and after 51.10 (non-dimensional) time from 
the beginning of the motion. All data necessary to implement equation (4.52) are 
available in [6]. 

The numerical tank has the following dimensions 


e tank length L = 20m; 
e water depth d = 1m; 


Some selected snapshots of the free surface evolution are presented in figures 4.19 
to 4.21. In particular, subplots of figure 4.19 show up to t = 25 (nondimensional 
time), while in figure 4.20 it is shown up to t = 49.5, which is just few instants 
before the steepness attains the critical value. Indeed, the three remaining subplots 
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Figure 4.16: Propagation of a regular wave generated by a piston wavemaker. 


shown in figure 4.21 depict the formation of the jet. From the figure it is also possible 
to note that the mesh has been refined a lot in the area where the cusp is forming. 
From this instant on it makes more sense to follow the evolution of the plunging 
breaker in an undistorted scale as in figure 4.22. 

To validate the code and its capability of correctly simulating fully nonlinear 
waves, a comparison with experimental results obtained in [6] is presented in the 
following. Six numerical gauges are set at nondimensional distances of a = 3.17, 
5.00, 6.17, 9.17, 10.83, 11.83 from the paddle and the time histories there numer- 
ically evaluated are compared with the experimental measurements. Figure 4.23 
shows the free surface elevation over time for each probe. 

The initial number of boundary elements of the free surface was NE = 120, 
then they self-adaptively increased whenever mesh refinement was needed. On the 
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(c) Free surface at t = 3s. 


Figure 4.17: Propagation of a regular wave generated by a piston wavemaker. 


other boundaries there were NE; = NE, = 3, NE, = 30 quadratic elements. The 
initial time-step adopted to integrate the boundary conditions was dt = 0.05s. 

For the purpose of this thesis it is crucial to describe the water particles velocity 
during overturning. Figure 4.24 gives an example on how accurately it is possible 
to investigate such velocities with the present code. 
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Figure 4.18: Propagation of a regular wave generated by a piston wavemaker. 
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Figure 4.19: Propagation of a wave packet generated by a piston wavemaker. 
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Figure 4.20: Propagation of a wave packet generated by a piston wavemaker. 
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Figure 4.21: Propagation of a wave packet generated by a piston wavemaker. 


63 


64 


03 


BEM Free Surface at t = 51.6445 s 


Hydrodynamic model 


BEM Free Surface at t = 51.6999 s 


03 


03 


0.1 


0.05 


(a) t = 51.64 


(b) t = 51.70 


ee Surface at t = 51.7499 s 


BEM Free Surface at t = 51.7999 s 


L 4 08h 
E 4 0.25} 
eer n= peer ET 

p AT »S | 0.2} pe 

_ N et 

hee NS Sew 
ee NO > set 
F 4 0.15 _ 
p 
yR k de 
ba a 
L SG: J oak 
Ma 
Sa 
— 
F ba 4 0.05} — 
h h f h h i J h h h h h 
15 11.6 117 11.9 12 115 11.6 117 11.8 11.9 12 
(c) t= 51.75 (d) t = 51.80 

BEM Free Surface at t = 51.8499 s BEM Free Surface at t = 51.8999 s 
r r 7 r 7 r T T r T r T 
L J oak 
F A 0.25} 

SE era ae a 
t a Ds | og% gee Sk 
A Pa 
_ SITÀ = N 

e: NS = "a / N 

P > pr 
F ( SN 4 cast a ( \ 

ui \\ a Ne N 
PG NV ae TT N 
Le I 4 O1p ST gi 
Tx I 

L ia 0.05} 
4 x 2 1 L L 4 ok L L L L L 
15 11.6 ILL 118 119 12 11.5 11.6 11.7 11.8 11.9 12 


(e) t = 51.85 


Figure 4.22: 


(f) t = 51.90 


Evolution of the plunging breaker. 
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Figure 4.23: Numerical and experimental time histories of the free surface elevation 


at six probes. 
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Figure 4.24: Velocities (red arrows) of the water particles at t = 51.70 of the spout 
evolution. Overturning wave generated by wave-wave interaction according to [6]. 
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4.5 Impact wave model 


The impact wave model adopted in the simulations performed in the next chapter 
is based on the Winke and Oumeraci’s work [8]. The model applies to overturning 
deep water breaking waves. Before introducing the impact model, a very short review 
of breaking waves is presented in the next section. After this, Morison’s equation 
will be recalled in order to prepare the background for the analysis of slapping wave 
loads. 


4.5.1 Breaking waves 


Among various special sea events present in nature, breaking waves definitively 
represent one of the most dangerous phenomenon. They may occur both in deep wa- 
ter and in shallow water, whenever the stability of the free surface is compromised 
by some reasons. In costal engineering the hazard represented by shallow water 
breaking waves is of primary interest, while when designing offshore structures deep 
water overturning breakers need to be taken into account carefully due their de- 
structive potentiality. The shallow water instability is reached when the ratio wave 
height /water depth attains the theoretical value of 0.78 [5]. 

While the deep water breaking condition is reached when the wave steepness 
en = H/L increases above the theoretical limit of 0.142. In the case of finite depth 
d, as it will be our case, the breaking limit becomes e, = 0.142 tanh (kd) where k 
is the wave number. When a water wave becomes unstable it may break in four 
different ways, giving rise to the so called spilling, plunging, collapsing, surging 
breaking wave profiles, as shown in figure 4.25. 

The most dangerous for the safety of offshore structures is definitely the plunging 
type which can cause strong and potentially destructive impacts. 

Currently, there is no systematic methodology to take into account such extreme 
events when designing offshore structures and, above all, unless using computa- 
tionally heavy CFD codes, they are only reconstructed by using empirical formula 
mostly suitable for deterministic design approaches. On this point it is worth quot- 
ing S.K. Chakrabarti who, in Chap. 3 of [10], says: 


“The theories described earlier for regular waves, including nonlinear 
Stokes waves and stream function theory, do not predict the kinematics 
and dynamic properties of very steep waves well. These waves are not 
only vertically unsymmetric, but also have large horizontal unsymmetry. 
If the design is based on these single steep waves, then a numerical theory 
need to be utilized. There are current attempts in describing such waves 
by the numerical wave tank methods and the method of New Waves [see 
Tromans, et al (1991) and Kim, et al (1999) for details]. These methods 
have not reached the design stage yet and are not commonly used in the 
design of offshore structures”. 


Plunging breakers definitely belong to the group the above quotation refers 
to when it mentions the horizontal asymmetry. Hence, what really deserves to be 
pointed out here is the fact that this thesis attempts to overcome what said in 
the last sentence of the above quotation. In fact, an efficient and computationally 
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Figure 4.25: Different types of breaking waves, source [5]. 


effective numerical tool able to simulate very steep waves in a systematic way, 


overturning plunging breakers included, is offered to improve the standard design 
procedures. 


4.5.2 Morison’s equation 


Wave forces on structures depend on many parameters, among which there are: 
time t, mean wave period T,, mean wave length L, horizontal water particles velocity 
v, member diameter D, kinematic viscosity u, member roughness x, etc. All of these 
variables can be rearranged and somehow normalized to define some conventional 


quantities 
e Keulegan—Carpenter’s number: KC = vT/D; 
e Reynolds’ number: Re = vD/y; 


e Diffraction parameter: D/L; 


e Relative surface roughness: e = «/D; 
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Morison’s equation applies only to those members which are considered slender, 
that is with the diffraction parameter lower than 0.2. Generally, we have 


D < 0.2 Morison’s equation 
L > 0.2 Diffraction/Radiation theory 


Note that the KC' number is also related to members main dimension. Consider 
in fact that for a fixed location, at the SWL, we have Nmax = H/2 and Umax = 
(gk/w) nmax- Recalling also the dispersion relation in deep water w? = kg, the KC 
number becomes 


H 
KC =r— 4. 
C TH (4.53) 


This means that the larger KC, the smaller the diffraction parameter. A general 
rule suggests that Morison’s equation is valid for KC > 6 [14]. 

To develop the impulsive wave impact model it is necessary to start by analyzing 
how Morison’s equation is composed [66]. The model takes into account the drag 
and inertia contributions induced by an undisturbed wave on a cylindrical mem- 
ber. Morison’s equation is a useful and simple engineering tool because it has the 
great advantage of considering the whole wave kinematics as if the cylinder were 
absent. On the other hand it presents some important limitations and drawbacks. 
For example it does not take into account wave run-up and it is questionable its 
validity for all wave theories, especially for highly nonlinear cases. In particular, by 
adopting Morison’s equation alone, extreme events like breaking waves cannot be 
considered at all. 

The force per unit-length given by Morison in [66] is the following 


f(t) = fo (t) + fm (t) = 
; pCp Dv (t) |v (t)| + i pC D754 (t) (4.54) 


where Cm and Cp are the mass and drag coefficients, respectively; py is the water 
density. 

The drag component in equation (4.54) would actually vanish if we had strictly 
to follow the the assumption of potential flow theory. Indeed, by virtue of this theory, 
one would end up into the famous D’Alebert’s paradox. Nevertheless, by assuming 
the real pressure distribution around a cylinder, once accounted for the dependency 
on Reynolds’ number, it is possible to describe the real drag force on the cylinder. 

Concerning the inertial contribution in equation (4.54), by virtue of the pure 
potential theory one would find Cm = 2. Looking at the inertial term fm we realize 
that (D/2)° x = vol is the unit-length volume of the cylinder, and thus 


fm (t) = 2p vol è (t) (4.55) 


which is nothing but twice the inertial force associated with a “slice” of the cylinder 
under consideration. In the reality it has been proved that the inertial coefficient 
is not really equal to 2 but it is given as Cm = 1+ km. The sub coefficient “1” 
represents the inertial force owned by the flow able to move an amount of fluid 
having the volume vol, which is actually replaced by the material member. This 
contribution assumes an undisturbed flow. 
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However, the presence of the cylinder modifies the flow forcing the water particles 
to encircle the pile. Thus, particles undergo an acceleration which is provided by a 
force called added mass. km represents the added mass coefficient and depends on 
the shape of the object. See [67] and [33] for further details. 

There are many researchers who proposed values for these hydrodynamic coef- 
ficients and of course it is not our intent to comment or list all of them. We just 
mention very few in order to get a practical idea. According to [13] typical values 
for the drag and inertial coefficients are given in table 4.2. 


| Re < 10° Re < 10° 
KC | <10 >10 <10 > 10 
Cp | 1.2 1.2 0.6 0.6 
Cm 2 1.5 2 1.5 


Table 4.2: Cp and Cw proposed in [13]. 


Values proposed by the American Petroleum Institute API and by the Society 
of Naval Architects and Marine Engineers (SNAME) are shown in table 4.3 


| smooth rough 
API SNAME API SNAME 
Cp | 0.65 0.65 1.05 1 
Cm | 1.6 2 1.2 1.8 


Table 4.3: Cp and Cm proposed by API and SANME. 


Coefficients recommended by Det Norske Veritas (DNV) in [7] are only related 
to the KC number and the relative surface roughness e. See figure 4.26. 

However, as most of the recommended values in this work follow instructions 
given in [2], reference values for hydrodynamic coefficients may be found also in the 
International Standard ISO 13819 - 2, Part 2: Fixed Steel Structures. 

The total force F = F (t) acting on the member is then obtained by integrating 
the unit length force along the member up to the free surface elevation as follows 


F(t) = Fp (t) + Fu (t) = 
ni Tage 
= J gPCpDu (t) lv (t) |dy + | > geCuD*é (t) dy (4.56) 
To perform the integral in equation (4.56) it is necessary to extend the solution 
of the wave kinematic model up to the free surface, therefore the commonly linear 
and weakly nonlinear theory need to be adjusted because they provide solution 


up to the SWL. In the numerical model developed in the next chapter, Wheeler 
stretching is usually adopted, [68] and [10]. 


4.5.3 Impulsive load due to plunging breakers 


Impact forces acting on offshore wind turbines can be two to four times larger 
than the non-impact forces stemming from waves of similar amplitude. And the 
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Figure 4.26: Hydrodynamic coefficients recommended by DNV, [7]. 


pressure due to impacts may be ten times larger than the non-impact pressure and 
it rises in a fraction of second. 

Traditionally, the contribution of slamming waves is accounted for in a deter- 
ministic sense without considering its temporal development. Indeed, the drag com- 
ponent of Morison’s equation (4.56) is just amplified by a factor which has a wide 
range of variation, typically from 2.5 to 5.15 [8], [69], [10]. [70]. 

Without considering the real time-history (though very short) of the impulsive 
contribution it is of course not possible to integrate such type of action in a full 
time domain analysis which accounts the global dynamic behaviour of the system. 


Contrarily, the slapping contribution due to plunging breakers can be described 
by adjusting the original Morison’s equation as follows 


F (t) = Fp (t) + Fu (t) +F (t) (4.57) 
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where F is the impulsive term. Wienke and Oumeraci in [71] and [8] developed 
a new analytical model to determine the impact force which results very accurate 
such as to represent a valid alternative to the earlier models proposed in [72], [32], 
[73], [74], among others. 

Adopting a potential flow model and neglecting the surface tension as well as 
the the forces due to gravity, they found that for a rigid cylinder of radius R the 
impulsive force for unit length in 2D is given by 


fr (t) = pRv* cosy [2r cos. — 2 cos tl . 


1 1 v 
. arctan 1- leon Di (4.58) 
for TE. 
O<t< -= — 4.59 
T ~ 8cosy v ( ) 
And 
T 
1 1 8 v a 
t = 2 t n 
fr (t) = pRv* cosy lion (Zos ) | 
arctanhy/1— Zy È Zy (4.60) 
R cosy R 
o 3 1R 12 1 R 
SES (4.61) 
32 cos y V 32 cos y V 
where i ay ae 
t=t- = 4.62 
32 cosy V ( ) 


The duration T; of the impact considered from the immersion time of the front 
line to the time of complete immersion of the half ellipse, see figure 4.27, is given 
by the following equation 

13 1 R 


t 32 cosy v 


(4.63) 


Figure 4.28 shows the fundamental sketch the above formulation is referred to. 
4.5.4 Numerical treatment of the plunging jet 


The subroutine implementing equations (4.58) and (4.60) will need, as passed 
variables, both v and m. These two values are obtained from the numerical simulator 
discussed earlier. To this end, the forming water jet has to be analyzed carefully. 

We start saying that for each time step of the simulation the free surface T ș (t) is 
numerically known. So that in the space-time neighborhood of the expected impact 
event, a dedicated subroutine finds at which time the impact would happen. Let us 
imagine that the turbine is located at x, = 35 m, see figure 4.29. 
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Figure 4.28: Sketch of the wave impact model. Image from [8]. 
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Figure 4.29: Example of the imminent overturning wave hitting the structure. At 


this time, n and v are computed. 


At this impact time, the free surface looks like the example shown in figure 4.29. 
To get the maximum wave height me is trivial, while to compute the impact velocity 
v an “averaging region should be identified”. 

After many numerical runs we discovered that averaging the horizontal velocity 
components over a region included by the maximum wave elevation node and the 
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node being the turning point for the free surface curvature, gives acceptable results 
in estimating the impact velocity v. The turning point is found internally by the 
subroutine by checking for each element when a change of the curvature sign occurs’. 


See figure 4.29. 


3For the j-th boundary elements the curvature k%) is computed by using the well known 
formula i 3/2 
k9) = [(ey — 92) /(2? + 9°)] 
where the derivatives are trivially carried out due to the simplicity of the shape functions. 


Chapter5 
Coupled wind—fully nonlinear waves model 


This chapter represents the second main part the thesis, where the fully nonlinear water 
waves simulator developed in Chapter 4 1s coupled with a hydro-aero-elastic solver of the 
whole system. The chapter begins by describing the solvers, then it enters into the details of 
the wind and wave loads simulations, and finally the fully coupled analyses are presented. 


5.1 Solvers description 


Here we will introduce the basic features of FAST: a combined modal and multi- 
body dynamics simulator [9]. The version used in this thesis is FAST_v602c-jmj 
which is a not-yet-released version. As it is an alpha version all the the new capa- 
bilities of the software, in particular the hydrodynamic module, are not yet docu- 
mented. 

However, due to the support provided by the National Renewable Energy Labo- 
ratory (NREL, Colorado)', it was possible not only to use the new modules but also 
to set up and implement the new impact model developed in this thesis. A detailed 
description of this new capability is presented in the next section. 

For three-bladed wind turbines FAST has up to 24 DOFs. All DOFs come from 
modeling rigid and flexible system components: tower, blades, etc. Reference [9] 
offers a detailed description of each DOF. The modus operandi of FAST (the parts 
relevant for our case) is sketched in figure 5.1, which also shows the input/output 
structure of the code. 

The aeroelastic forces acting on the rotating blades are computed by AeroDyn 
which is internally called by FAST’s main program at each time step. To compute 
the aerodynamic forces, AeroDyn implements both the Blade Element Momentum 
theory and the Generalized Dynamic Wake model [35], see Chapter 3. 

To give an overview of the fundamental features of the solver we will shortly 
describe all the input files which are necessary to the simulations, see figure 5.1. 


5.1.1 FAST input 


e Primary input file: *.fst: contains many sets of the parameters which have 
to be entered. We mention just the following: 


— PLATFORM which calls an additional input file containing the platform prop- 
erties (.dat) 


In the person of Jason Jonkman. 


Enzo Marino, An integrated nonlinear wind-waves model for offshore wind turbines 
ISBN 978-88-6655-051-8 (print) ISBN 978-88-6655-053-2 (online) © 2011 Firenze University Press 
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Platform 


Time series 


Figure 5.1: Selected files involved in FAST and AeroDyn. 


TOWER which calls an additional input file containing the tower properties 
(.dat) 

— BLADE which calls additional input file(s) containing data about the blades 
(.dat) 

AERODYN which calls an additional input file containing AeroDyn input pa- 
rameters (.ipt) 


* the .ipt file calls the file containing the wind data (.wnd) 
* the .ipt file calls the file containing the airfoil data (.dat) 


OUTPUT which itemizes the output channels required 


e Platform (monopile) input file *.dat: also in this case we mention the 
blocks more relevant for our scope: 
— MASS AND INERTIA (relative to the platform) 
— PLATFORM (loading model 0: none, 1: user-defined from routine UserPtfmLd 
(switch)) 
TOWER (Loading model) 
— WAVES 
— CURRENT 
OUTPUT 


e Tower input file *.dat: the data sets required in this file are: 
— TOWER PARAMETERS 
— TOWER ADJUSTMUNT FACTORS 

TOWER FORE-AFT MODE SHAPES 

TOWER SIDE-TO-SIDE MODE SHAPES 


e Blade input file *.dat: 
— BLADE PARAMETERS 
— BLADE ADJUSTMENT FACTORS 
— DISTRIBUTED BLADE PROPERTIES 
— BLADE MODE SHAPES 
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5.2 5-MW Baseline reference model 


This section is devoted to a brief description of the model adopted in all the 
simulations carried out to test the new hydrodynamic model proposed. 

The turbine model is the 5-MW Reference Wind Turbine for Offshore System 
Development. See Figure 5.2 for its general layout. 


wind 


126m 


tower 90m 


SAWA mono pile ‘Lg 
WS support | 


Figure 5.2: Layout of the model. 


This model has been created by J. Jonkman et al. at the National Renewable 
Energy Laboratory (NREL, Colorado) with the aim of creating a common reference 
basis to pave the way for further investigations concerning offshore wind turbines lo- 
cated both in shallow and deep waters. In [75] the authors specified all the technical 
characteristics of the hypothetical, but realistic, multi-megawatt large wind turbine 
model by combining some available data from different machines (e.g. REpower 5- 
MW and Arveva Multibrid 5000) together with data assumed in some conceptual 
models of projects like RECOFF, etc. NREL 5-MW Baseline Wind Turbine gross 
properties are itemized in table 5.1. 

In the following sections a short description of the model is given. Further details 
are available in [75], [76] and their bibliography. 


5.2.1 Rotor and support structure 
Rotor blades 

The turbine has three blades with structural properties of the 62.6 m—long LM 
Glasfiber blade used in the DOWEC study. As this type of blade is 1.1m longer 


than the 61.5 m-long LM Glasfiber blades adopted on the actual REpower 5-MW 
machine, in [75] they truncated the 62.6 m—long blades at 61.5 m span to obtain the 
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Rating Power 


5 MW 


Rotor Orientation, Configuration 


Upwind, 3 Blades 


Control Variable Speed, Collective Pitch 
Drivetrain High Speed, Multiple-Stage Gearbox 
Rotor, Hub Diameter 126m, 3m 

Hub Height 90m 


Cut-In, Rated, Cut-Out Wind Speed 


3m/s, 11.4m/s, 25 m/s 


Cut-In, Rated Rotor Speed 


6.9rpm, 12.1rpm 


Rated Tip Speed 80 m/s 
Overhang, Shaft Tilt, Precone 5m, 5°, 2.5° 
Rotor Mass 110t 
Nacelle Mass 240t 

Tower Mass 347.460t 


Coordinate Location of Overall CM 


(—0.2 m, 0.0 m, 64.0 m) 


Table 5.1: Key properties of the NREL 5-MW Baseline Wind Turbine. 


structural properties of the NREL 5-MW baseline blades. Then properties at 61.5 m 
have been found interpolating the two properties at 61.2 m and 61.7 m stations. More 
exhaustive data are available in [77]. Each single blade is made up of three different 
varying geometries. Each different geometry is described in its special file where the 
aerodynamical properties are assigned. See table 5.2. Appendix B of |75] reports all 
the mentioned files. 


Tower and pile 


The tower adopted in the NREL baseline model is tubular shaped and has 
the same geometric characteristics of the tower used in the DOWEC study [76]. 
Table 5.3 summarizes the key geometric parameters. The tower diameter and wall 
thickness are assumed to vary linearly. The base diameter of 6 m is equal to the 
diameter of the monopile. 

Regarding the substructure, the main properties are itemized in table 5.4. 

Note that in general FAST uses the following variables to model the platform. 
As shown in Figure 5.3, TwrDraft denotes the downward distance from the Still 
Water Level (SWL) to the tower base platform connection. PtfmCM and PtfmRef 
represent the distance from the SWL and the platform center of mass and the 
platform reference point, respectively. The latter is the point where the platform 
DOFs are located. These three parameters are useful for modeling floating platforms, 
in our case they all equal the water depth d = 20 m. 
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Figure 5.3: Platform scheme. Image from from [9] 
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Node RNode AeroTwist ARNode Chord Airfoil Table 

1 2.8667 13.308 2.7333 3.542 Cylinder1.dat 

2 5.6000 13.308 2.7333 3.854 Cylinder1.dat 

3 8.3333 13.308 2.7333 4.167 Cylinder2.dat 

4 11.7500 13.308 4.1000 4.557 DU40_ Al7.dat 

5 15.8500 11.480 4.1000 4.652 DU35_A17.dat 

6 19.9500 10.162 4.1000 4.458 DU35_A17.dat 

T 24.0500 9.011 4.1000 4.249 DU30 _A17.dat 

8 28.1500 7.795 4.1000 4.007 DU25_ Al7.dat 

9 32.2500 6.544 4.1000 3.748 DU25 _A17.dat 

10 36.3500 5.361 4.1000 3.502 DU21 Al7.dat 

11 40.4500 4.188 4.1000 3.256 DU21 Al7.dat 

12 44.5500 3.125 4.1000 3.010 NACA64 A17.dat 
13 48.6500 2.319 4.1000 2.764 NACA64 A17.dat 
14 52.7500 1.526 4.1000 2.518 NACA64 A17.dat 
15 56.1667 0.863 2.7333 2.313 NACA64 A17.dat 
16 58.9000 0.370 2.7333 2.086 NACA64 A17.dat 
17 61.6333 0.106 2.7333 1.419 NACA64 A17.dat 


Table 5.2: Distributed blade aerodynamic properties. 


Tower base diameter, wall thickness 6m, 0.027 m 


Tower top diameter, wall thickness 3.87 m, 0.019m 


Table 5.3: Tower geometric properties. 


Pile length?, diameter 20m, 6m 


Pile wall thickness, total weight 0.060 m, 187.90t 


Table 5.4: Monopile properties. 


5.3 New slamming wave module in FAST 


As already mentioned, FAST is a Fortran open source code developed at NREL 
with the extraordinary feature of having some subroutines which can be defined by 
the user upon their specific needs. This was, indeed, the case of this work where the 
impact forces associated with some probable extreme seas, first computed externally 
by dedicated simulations, are passed to the time marching solver by means of the 
used defined subroutine UserTwrLd. This subroutine provides user defined tower 
loading in the case of monopile substructure. 

Before describing how the impact model has been implemented, it is useful to 
give just an overview about the subroutine MorisonTwrLoading. The latter is stored 
in the HydroCalc.f90 file and provides the forces acting on the tower calculated 
by means of Morison’s equation. File Hydro Calc.f90 also contains all the necessary 
subroutines to generate the requested sea: linear regular or linear irregular. The 
latter can be generated by adopting either the Pierson—Moskowitz or the JONSWAP 
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spectrum. However we will always use our external sea generator. 

The most important input parameters that control the wave kinematics are pro- 
vided by the “platform monopile” input file (see previous section). In the present case 
the following parameters are frequently invoked: TwrLdMod (from block TOWER), 
which sets the tower loading model (0: none, 1: Morison’s equation, 2: user-defined 
from routine UserTwrLd). 

From block WAVES, the parameters WaveMod, WaveStMod WaveT Max, WaveDT 
are crucial to develop the new impact model. The first, WaveMod, sets the type of 
incident wave kinematics (0: none, that is still water, 1: plane regular, 2: irregu- 
lar with JONSWAP or Pierson—Moskowitz spectra, 3: user-defined spectrum from 
routine UserWaveSpctrm, 4: GH Bladed wave data). WaveStMod sets the type of 
stretching, that is how to extend the linear solution, which is valid only up to 
the still water level, up to the instantaneous free surface elevation. Switches for 
WaveStMod are 0: no stretching, 1: vertical stretching, 2: extrapolation stretching, 
3: Wheeler stretching. WaveTMax and WaveDT set the total simulation time and 
the time step, respectively. 

To implement the new slapping wave module we first need to pass to FAST the 
wave kinematics consistent with the actual extreme sea which may contain some 
breaking wave events. Hence, the random sea is generated externally, see section 5.4 
for details, then it is passed to FAST by the same way it reads in the wave kinematics 
files generated by GH Bladed. Namely, three files are necessary. They all have the 
same name root which can be assigned as input. 


e FNL_ FAST. tat 
e FNL surface.tat 
e FNL_kinematics.tat 


where FNL (in this case) is the passed name root. 

FNL_FAST.tat contains the coordinates of the points where the kinematics is 
provided. The coordinates are expressed with respect to a planar system with the 
vertical axis y upward oriented and passing by the center of the monopile, while 
the x-axis points in the direction of the wave propagation. The origin is fixed 
at the still water level. FNL_ surface.trt stores the free surface elevation at the 
wind turbine location (x; = 0) during the whole simulation time WaveT max, while 
FNL_kinematics.tat provides: 


e particles velocity along x-axis , v” 

e particles velocity along z-axis , v (always zero in 2D model) 
e particles velocity along y-axis , v” 

e particles acceleration along x-axis , a” 

e particles acceleration along z-axis , a? (always zero 2D model) 


e particles acceleration along y-axis , a” 


e dynamic pressure pp 
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To understand the structure that FNL_ kinematics.tat must have to be correctly 
read in by FAST, here there is an example of the Matlab function generating it: 


function [] = write_kinematics(Elevx0,time,y) 


global Aw omega k_w phase rho_w g d 

% NOTE: xt = 0! 

fid = fopen(’fast\FNL_kinematics.txt?,°w?); for i = 1:length(time) 

for h = 1:length(y) 

% With Wheeler stetching 
eta = Elevx0(i); 
if y(h) <= eta 
q = d/(dteta); 
vx = sum(Aw.*omega .* cosh( k_w.*( d + q*y(h) + d*(q-1)) )./sinh(k_w.*d).*cos(-omega.*time(i)+phase)); 
vy = sum(Aw.*omega .* sinh( k_w.*( d + q*y(h) + d*(q-1)) )./sinh(k_w.*d) .*sin(-omega.*time(i)+phase)) ; 
ax = sum(Aw.*omega.72 .* cosh( k_w.*( d + q*y(h) + d*(q-1)) )./sinh(k_w.*d).*sin(-omega.*time(i)+phase)); 
ay =-sum(Aw.*omega.72 .* sinh( k_w.*( d + q*y(h) + d*(q-1)) )./sinh(k_w.*d).*cos(-omega.*time(i)+phase)); 
p = rho_w*g*sum(Aw .* cosh( k_w.*( d + q*y(h) + d*(q-1)) )./cosh(k_w.*d) .*cos(-omega.*time(i)+phase)) ; 
else 


fprintf (fid,?%d Yad Yad Yad Yad Yad %a\r\n?,vx, 0, vy, ax, 0, ay, P); 
end 
end fclose(fid); 


To check whether the wave kinematics provided by using external text files were 
correct, a very simple test comparing results obtained adopting both simulators has 
been carried out. 

A simple regular wave has been generated first by using the internal wave solver, 
then by passing data provided by our external solver. For the sake of simplicity no 
wind is assumed to blow and the turbine is parked. The incident regular periodic 
wave is characterized by the following wave height and wave period: H = 12m and 
T = 12s. 

For both the simulations shown in figures 5.4 and 5.5, hydrodynamic forces 
are computed by using Morison’s equation, i.e. TwrLdMod: 2, while WaveMod is 1 
when the wave is simulated internally and 4 when read in from text files. Note that 
WaveMod 4 means loading GH Bladed wave files, but in the present case our own 
simulator is used. Only the file format and their organization is passed according 
to GH Bladed file type. 

By comparing the tower base forces shown in figure 5.4(a) and 5.5(b), it is 
evident that the kinematics passed gives the same results than WaveMod 1. In fact, 
forces have the same period and amplitude. 

Moreover, as expected, the same agreement is observed by comparing the tower 
base moments. Indeed, also for them, the structural response shown in figures 5.5(a) 
and 5.4(b) is basically the same. 

Note that both for tower base forces and moments, meaningful comparisons are 
made between Mys and Fyt. The remaining internal forces are reported just for the 
sake of completeness. 

At this stage, once it has been proved that the external wave solver developed in 
the framework of this thesis works well, it is possible to move on to the impact model. 
As already said, a short introduction of MorisonTwrLd may be helpful as Morison’s 
equation is the starting point for the impulsive contribution, see section 4.5. The 
subroutine needs as input variables the current tower node, the tower diameter, 
the inertial and drag coefficients, the three components of the translational and 
the three components of the rotational displacements, the three components of 
the translational and the three components of the rotational (angular) velocities. 
Moreover, the current simulation time is also needed. 
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Figure 5.4: Tower base forces obtained with wave kinematics computed internally 
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(a) Tower base forces computed by using FAST’s internal wave 
simulator. 
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(b) Tower base forces computed by using the external wave kine- 
matic solver and passed via text files. 


by FAST (WaveMod: 1) and passed form outside (WaveMod: 4). 
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MorisonTwrLd subroutine outputs the force vector TwrFt, which is formed by the 
x (surge), y (sway) and z (heave) components of the portion of the tower force per 
unit length at the current tower element and the roll, pitch, and yaw components 
of the portion of the tower moment per unit length acting at the current tower 


element. 
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(a) Tower base moment computed by using FAST’s internal 
wave simulator. 
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(b) Tower base moment computed by using the external wave 
kinematic solver and passed via text files. 


Figure 5.5: Tower base moments obtained with wave kinematics computed internally 
by FAST (WaveMod: 1) and passed form outside (WaveMod: 4). 


5.3.1 Slamming tower loading subroutine 


The slamming tower load subroutine developed in this thesis adds a new feature 
to FAST. In fact, whenever under extreme climates a designer detects that highly 
nonlinear wave are present in the simulated sea and then the risk of impacts is pretty 
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hight, it is possible to switch the input parameter TwrLdMod to 2 in order to call 
the special user-defined subroutine implementing the theoretical model discussed 
in section 4.5. 

It was necessary to introduce MorisonTwrLd subroutine because our user defined 
one contains a call to it. This should not surprise as the impulsive contribution is 
just added to the drag and inertial ones already provided by Morison’s equation [8]. 

So the idea is: if breaking wave events are expected, then switch TwrLdMod to 
2, else switch TwrLdMod to 1. Algorithm 2 summarizes this logic. 


Algorithm 2: Basic logical scheme on how the variable TwrLdMod makes 
FAST compute tower loads. 

input : TwrLdMod 

output: TwrFt 


if TwrLdMod = 1 then 
No impacts are expected, TwrFt = fp (t) + fm (6) 
| MorisonTwrLd is called; 
else if TwrLdMod = 2 then 
Impacts may occur, TwrFt = fp (t) + fm (© +fr(1) 
| UserTwrLd is called; 


According to algorithm 2, when UserTwrLd is called by switching TwrLdMod to 
2, the first thing the routine does is to read in the following text files: 


o NtNb.tat 
e eta_b.tat 
© curl.tat 
o Tb.tat 

e f_ Impact 


To understand the logic of the routine, the above files need to be shortly de- 
scribed. NtNb.tat just contains the two parameters: Nt and Nb, which denote the 
number of total time steps concerning the sea simulation (Nt = Tsim/WaveDT + 1) 
and the number of expected impact events in that sea state, respectively. 

File eta_ b.txt provides the maximum wave heights during the impacts and of 
course the file consists in Nb values. To remind the importance of the variable my 
(known as eta_b in the code) it is useful to recall the impact scheme proposed by 
Wienke and Oumerci in [8] and shown again in figure 5.6. 

To properly compute the impact load, it is also extremely important to define 
the so called impact area. See figure 5.6. This can be done by means of the curling 
factor A which basically gives the portion of tower being hit by the overturning 
spout of water. It would be very interesting in future to use the fully nonlinear 
solver here developed to asses and compare with experiments this factor, but in this 
application it has been fixed according to literature as A = 0.46 [78], [79], [8]. 
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area of 
impact 


Figure 5.6: Sketch of the wave impact model. Image from [8]. 


File Tb.txt provides the vector gathering the times at which the impacts are 
expected. During an extreme sea state, in fact, several plunging breakers may occur 
at the structure location, thus, to take into account all of them, this variable needs 
to be known?. 

Finally, reading in the file f_ Impact, the subroutine will create an allocatable 
variable called f_Impact. This variable has Nt rows and three columns: time, force 
in x direction (main direction of the impact), force in y direction. Of course the 
third column is always zero for the present model. 

Now, when all the input data have been read in by the routine, the fraction of 
the current tower element included in the impact area is computed (this is necessary 
because the impact force applies of course only to those elements, or portion of them, 
included in the range Am). To do this a logic similar to that used in MorisonTwrLd to 
find out which elements fall between the mudline and the instantaneous free surface 


has been implemented. 
A selection of the most relevant lines of UserTwrLd are listed in the following: 


!!!!Some key lines of: SUBROUTINE UserTurLd ( JNode, X, XD, ZTime, 

DirRoot, TwrAM, TwrFt ) ! Select eta_b and curl corresponding to the 

current impact event 

CALL MorisonTwrLd ( JNode, TwrDiam, TwrCA, TwrCD, X, XD, ZTime, TwrAM, TwrFt ) 


! Initialize eta_b and curl 


eta_b 
curl 


0 
0 


DO I = 1,Nb ! Loop through the impact events 
IF ( ( Ztime >= Tb(I) - 1) .AND. ( Ztime <= Tb(I) + 1 ) ) THEN 


eta_b = eta_b_vec(I) 
curl = curl_vec(I) 
END IF 


END DO 


3Note that this file provides the variable Tb to the subroutine, and it is nothing but what we 
will call # later on in section 5.5 of this chapter. 
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IF ( ( WaveKinzi0(JNode) - 0.5*DZNodes(JNode) ) >= eta_b ) THEN 
DZFractW = 0.0 
ELSEIF ( ( WaveKinzi0(JNode) + 0.5*DZNodes(JNode) ) <= eta_b ) THEN 
DZFractW = 1.0 
ELSE 
DZFractW = ( ( eta_b - ( WaveKinzi0(JNode) - 
0.5*DZNodes(JNode) ) )/DZNodes(JNode) ) 


ENDIF 

IF ( ( WaveKinzi0(JNode) - 0.5*DZNodes(JNode) ) >= (1.0-curl)*eta_b ) THEN 
DZFractS = 1.0 

ELSEIF ( ( WaveKinzi0(JNode) + 0.5*DZNodes(JNode) ) <= (1.0-curl)*eta_b ) THEN 
DZFractS = 0.0 

ELSE 
DZFractS = ( ( ( WaveKinzi0(JNode) + 0.5*DZNodes(JNode) ) - 

( (1.0-curl)*eta_b ) )/DZNodes(JNode) ) 
ENDIF 


DZFract = DZFractW*DZFractS 


! Compute the impact load (for the moment it is simply read from 
outside): 


IF ( DZFract > 0.0 ) THEN 


f_imp(1) = InterpStp ( ZTime, f_Impact(1,:), f_Impact(2,:), LastInd, Nt) 
f_imp(2) = InterpStp ( ZTime, f_Impact(1,:), f_Impact(3,:), LastInd, Nt) 
DO K = 1,2 ! Loop through the xi- (1) and yi- (2) directions 


TwrFt(K) = TwrFt(K) + f_imp(K)*DZFract 
ENDDO ! K - The xi- (1) and yi- (2) directions 


ENDIF 


It is important to observe that in the subroutine UserTwrLd, first MorisonTwrLd 
provides the drag and inertial terms, then the remaining part adds the impulsive 
contribution. In addition to this, two important things should be pointed out from 
the above subroutine: first, it makes use of an additional subroutine named InterpStp 
which provides the value of f impact interpolated at Ztime. This interpolation per- 
mits to compute the impact force with its own time step, which normally is WaveDT 
and then to be included in FAST main solver which has a different time step DT. 
It should also be considered that dealing with impulsive events, let us say with an 
average duration of magnitude 0.01s, to capture the impulsive load it is necessary 
that a minimum number of time steps have to be included in the impact duration 
time. This induces some restriction both on WaveDT and DT. On the contrary, to 
generate the wave kinematics a time step of 0.5s would be enough, but it could not 
capture the impact events. For this reason the impact force is always given with a 
proper time step. 


The last remark about the above routine (see the last lines) is that the interpo- 
lated values of f_Impact, f_imp(1) and f_imp(2) are added to Morison’s contribu- 
tions at each time step. This can be done because f_Impact is nonzero only for the 
duration of the impact(s). 
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5.3.2 Test. of the slamming load subroutine 


At this point of the work, all the necessary tools to set up a preliminary test 
of a slamming wave have been provided. Plots showed in figures 5.4 and 5.5 were 
the time histories of the structural response subjected to a weakly nonlinear regular 
wave characterized by a wave height H = 12m and a wave period T = 12s. This 
regular wave was rather far from the breaking limit, indeed its steepness is e = 0.079. 

In this section, on the contrary, we need to increase the wave steepness in order 
to generate a plunging breaker from which, by following the scheme described in 
section 4.5.4, the input data to be passed to subroutine UserTwrLd can be derived. 

A satisfactory high value of the steepness is achieved by doubling the wave height 
of the previous case, so that H = 24m, while the wave period is unchanged: T = 12s. 
Therefore the steepness becomes e = H/L = 24/152.365 = 0.1575. The free surface 
evolution with the formation of the overturning wave is shown in figure 5.7. 

With a tower diameter of 6m, the resulting impact force due to the breaker of 
figure 5.7 is plotted in figure 5.8. 

The impact force vector, after being written in f_ Impact.tat, is passed as variable 
f_ Impact to FAST and it has the classical impulsive form as shown in figure 5.10. 

According to this time history, indeed, the second column of file f_ Impact.tat 
is nonzero only in a proper neighborhood of ty = 12s. The scale of figure 5.10 is 
too large to capture the real distribution of the impulsive force (which is shown in 
figure 5.8) but, as already said, FAST integration time step has been chosen in such 
a way to sample the impulsive contribution in an enough number of points. 

Finally, it is possible to move to the structural response. As did earlier when we 
tested the wave kinematic solver, here we want to show the effect of the tower base 
internal forces, basically the shear force F,; (that is in the same direction of the 
incident wave) and the bending moment Myt rotating around the y-axis. The other 
internal forces and moments are also presented just to check that nothing occurs in 
the remaining directions. 

The first subplot of figure 5.11(a) clearly describes that the periodical response 
due only to one single-harmonic loading action is suddenly shocked by the impul- 
sive action associated with the slamming wave. The shear forces in the direction 
of the wave motion, indeed, presents a peak right at t = 12s, compare with fig- 
ure 5.10, while the force orthogonal to the main direction of the wave Fy; is not 
affected at all by this contribution. Another clear effect is registered by the time 
history of the vertical force F,; at the tower base. In this direction, see the lower 
subplot of figure 5.11(a), the impulsive force also gives a relevant effect because the 
structure, after being hit, starts oscillating and this activates the the rotor mass, 
which, together with the mass of the tower, starts exciting the structure also in 
the vertical direction giving rise to a transient response. After a few seconds such a 
transient behaviour is damped. The vertical force comes back to the structure self 
weight after approximately 25s, while a shorter damping period is necessary for Fyt 
to recover its periodicity. However, more details about the effects of the slam upon 
the vertical tower base force will be provided in the last paragraph of this section. 

Observing in depth F,;, however, a short delay between the peak due to the 
impact and the maximum fore due to Morison’s contributions alone occurs. This 
is due to the fact that the impact force is computed with a nonlinear solver, while 
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Figure 5.7: Free surface evolution of a steep regular breaking wave. Red arrows 
denote the free surface particles velocity and the blue dots the boundary element 


mesh. Input data from table 5.12. 
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Figure 5.8: Impact force per unit length associated 
in figure 5.7(c). The impact duration is T; = 0.094s. 


0.1 


with the breaking wave shown 


-10 0 10 


Figure 5.9: Closer view of the forming plunging breaker shown at the same time of 
figure 5.7(c). From this configuration the impact velocity and maximum wave height 
have been extracted in order tho get the time history of the impact load shown in 


figure 5.8. 
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Figure 5.10: Impact force time history stemming from a steep regular wave breaking 
at t= 12s, with m = 11.40m, A = 0.46. Total simulation time Tsim is 180s. 


Morison’s force stems from a linear wave model. This delay is probably even more 
increased due to an excessive use of smoothing, which tends to dissipate some energy. 
A more careful use of both smoothing and regridding subroutines has been done in 
the next simulations. 

Likewise, looking at figure 5.11(b), the results are pretty encouraging as neither 
the yaw moment nor the bending moment around the main direction of the incoming 
wave are affected by the hit, while, on the contrary, maximum effects is produced 
on the bending overturning moment. Here, in the middle subplot of figure 5.11(b), 
a clear peak occurs at 12s which dissipates approximately after 10s. Note that the 
damping period gives an idea of how much the whole structure suffers the shock. 

It is also interesting to observe the tower top displacement in the incident wave 
direction. This displacement is shown in figure 5.12. 

Also for the tower top displacement at t = 12s, when the impact occurs, an 
increase of the tower deflection is registered. The series plotted in figure 5.12 is too 
short to really compare the transient effect with the stationary behavior. It seems 
that after a peak of about 26cm due the slam, a progressive decay of the following 
peaks is registered. 


Further considerations on F, time-history 


To be sure that the above reasoning to justify the effect on F, is correct, the 
same simulation has been run disabling the tower degrees of freedom. Therefore, 
by assuming that the tower cannot undergo any displacement, that is no fore- 
aft oscillation is permitted, then no excitation is caused in the vertical direction. 
Evidences of this are given in figure 5.13. 

The explanation given above about the transient effect also in F,; time history 
is absolutely confirmed. Due to the absence of tower oscillations the impact effect 
is even clearer than before. The lower subplot of figure 5.13(a) shows that in this 
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(a) Tower base forces. 
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(b) Tower base moments. 


Figure 5.11: First test on the structural response accounting for the impulsive load 


generated by a plunging breaker obtained from a very steep regular wave. (WaveMod: 
4, TwrLdMod 2). 
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Figure 5.12: Tower top fore-aft displacement of the turbine parked and subjected 
only to hydrodynamic loads including the impact force associated with a plunging 
breaker. 


case no transient effect is induced in z-direction and, as expected, the tower base 
vertical force remains constant on 8578 kN (as after the transient of figure 5.11(a)) 
which represents the self weight. 

The last remark before closing this important test of the new impact model 
implemented in FAST regards the peaks of F,; and My: plotted in the upper and 
middle subplots of figures 5.13(a) and 5.13(b), respectively. 

The peak value of the shear force is Fyt,, = 1.8 x 10° kN, while the maximum 


shear force during the stationary phase is Fyt = 5.8 x 10° kN. Thus, the impact 
increases in the shear force by 210%. While the peak value of the bending moment, 
is Myt,, = 4.8 x 10 kNm and the maximum value during the stationary phase is 


My = 8.7 x 10° kNm. This means that due to the impact, the fore-aft bending 
moment undergoes an augmentation of 451%, that is Myt, = 5.5Myt. Note that 
these amplifications should be treated carefully at this stage of the work, because 
the plunging wave has been generated without taking into account a realistic sea 
state. Indeed, the aim of this section is only to test the efficiency of the numerical 
model. 
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Time Series Plots of Tower base forces 
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(b) Tower base moments. 


Figure 5.13: The same impact load of figure 5.8, but in order to better investigate 
the nature of Fzt, this figure shows the tower base forces and moments when the 
tower degrees of freedom are disabled. 
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5.4 Wind and wave loads generation 


According to TEC61400-1 [36], an Extreme Wind Speed Model (EMW) is as- 
sumed to estimate the wind loading conditions experienced by the offshore wind 
turbine. Consequently, as we are here mainly interested in extreme events, the rele- 
vant sea conditions are basically made up of wind-generated waves. So that, based 
upon the assumed wind, the most probable sea state around the offshore monopile, 
defined in terms of significant wave height and mean wave period, is reproduced in 
the time domain through a linear spectral approach. In particular the JONSWAP 
spectrum is adopted as the most suitable in case of wind—waves. For a generic time 
instant t the extreme scenario looks like the sketch in figure 5.14. 
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Figure 5.14: Sketch of the wind turbine located at x; = 0 in the 2D spatial domain 
Di = [tmin, Tmax] for a given time instant t. 


Figure 5.14 describes the main load sources for the whole structure acting to- 
gether and in a continuous mutual interaction. 

Once the wind-generated irregular sea is known in space and time, a suitable 
breaking wave criterion is used to check whether plunging breakers occur. And 
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thus, only when the steepness becomes large enough to induce wave instabilities, 
the fully nonlinear simulator is called in order to follow the sea evolution with an 
higher resolution in a narrower neighborhood of the substructure. This procedure 
is analyzed in depth on section 5.4.3. 

The wind speeds at each blade element at the current time are computed from 
the hub-height turbulent value by accounting only for the assumed wind shear. As 
the rotor lies in the parked configuration (all the three blades are set with 90° pitch 
angle), the conservative assumption of not considering a full field turbulent wind 
generated by means of coherence functions can be accepted. 


5.4.1 Wind loads 


As already mentioned an extreme turbulent wind speed model (EWM) is adopted 
according to IEC61400-1. The turbulent wind is generated at the hub—height by us- 
ing TurbSim which is an utility developed at NREL capable of generating both 
hub-height and full field wind data according to several spectral formulation. For 
our simulations the Kaimal model will always be adopted, see section 3.4 for details. 


Extreme wind speed generation 


In FAST’s primary input file (i.e. *.fst), in the AERODYN block, AeroDyn input 
file (i.e. *ipt) must be specified. In *.ipt (line 10) at least two additional input 
files have to be provided: the wind file, containing either the hub-height wind speed 
(steady or time varying) or a full-field data file. 

AeroDyn automatically detects the wind data type: if no file extension is specified 
(but only the root), then AeroDyn treats it like a full-field file and it expects to find 
both the files *.wnd and *.sum. Otherwise, a hub-height wind file, *.hh, is assumed. 
This will always be our case. To summarize the input data necessary to generate a 
turbulent hub-height wind file, a part of TurbSim input file is shortly commented 
in the following. 

TurbSim’s input file is made up of five blocks: 


e Runtime Options; 

e Turbine Model Specifications; 

e Meteorological Boundary Conditions; 

e Non-IEC Meteorological Boundary Conditions; 
e Coherent Turbulence Scaling Parameters; 


Parameters provided in the above blocks are related both to the turbine type 
being simulated and the Design Load Case requested. In particular, in the first block 
settings about pseudorandom number generation and output type are specified. 
In the second block, parameters concerning the grid size and density, hub height, 
analysis time, etc. are to be specified. Most of these parameters depend on the 
turbine type for which the wind is simulated. 

To set up the Kaimal model described in section 3.4, it might be opportune to 
show the forth block about “Meteorological Boundary Conditions”: 
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spbsssse Meteorological Boundary Conditions ------------------- 
"TECKAI" TurbModel\\ 

"1-ED3" IECstandard\\ 

acu IECturbc\\ 

"3WM50" IEC_WindType\\ 

default ETMc\\ 


"pL" WindProfileType\\ 
90 RefHt\\ 


37.5 URef\\ 


Even though offshore wind turbines fall in the wind class “S”, whose wind pa- 
rameters are specified by the designer (see Annex A of [2]), here for the sake of 
simplicity we assume the wind turbine is designed according with parameters valid 
for wind class III. Therefore we have U,e; = 37.5 m/s and for the turbulence inten- 
sity we choose the group “C”, that is I-er = 0.12, see [36]. As already mentioned 
in section 3.4, the simulations performed to test the coupled wind-fully nonlinear 
waves here developed, employ an extreme turbulent wind model. As in this case 
the standard deviation of the turbulence is always assumed to be depending only 
on the hub height wind speed, the turbulence intensity is not actually used in the 
following simulations. 

Given the hub height of 90m measured from the still water level, for the tur- 
bulent extreme wind speed model, the 10-min average wind speed profile with re- 
currence periods of 50 years and 1 year, respectively, are those given in section 3.4. 
The extreme mean wind profiles both for a 50 years and one year return period are 
shown in figure 5.15. 

Of course, such wind velocities are far beyond the cut-out wind speed, which is 
25 m/s, indeed in these two cases the turbine is parked or standstill. 

Note that in FAST, in order to set the turbine in a parked condition, it is 
necessary (i) to turn off all controllers; (ii) to put the blade pitch angles at 90° 
(basically no lift forces); (iii) deactivate the induction model. 


5.4.2 Wind-correlated sea states 


An extreme sea state is defined by its significant wave height and mean zero- 
upcrossing wave period which should be estimated by taking into account the actual 
correlation between the whole environmental processes. As already discussed in 
Chapter 2, this requires a joint probabilistic model for the weather parameters 
relevant to the problem under consideration: wind speed, significant wave height, 
mean wave period. 

As described in Chapter 2, the most suitable tool to this end is the so called 
Inverse FORM (First Order Reliability Method). Such a tool provides an envi- 
ronmental contour defining, for a required return period, environmental variables 
combinations at which the most extreme response lies. In other words, it is necessary 
to search along the contour in order to determine the point at which the conditional 
expected extreme response becomes the most extreme. The idea of investigating a 
contour, hence all combinations, is particularly effective because in most cases the 
critical environment is structure-dependent. For instance, assume to have found a 
contour of significant wave height and current velocity, then in the case of shallow 
water, it will be more dangerous the combination which maximize the significant 
wave height rather than the one which maximize the current velocity, [80]. 
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EWM profile, U,.; = 37.50 m/s 


0 10 20 30 40 


Figure 5.15: Extreme Wind Model velocity profiles. 


A joint probabilistic model of mean wind speed, significant wave height and 
spectral peak period needs to be built. The wind speed is chosen as the primary pa- 
rameter, the significant wave height is assumed to have second most influence. The 
joint model is used to establish a contour surface, giving combinations of the envi- 
ronmental variables with return period of 50 years according the recommendations 
in [2]. 

The reliability of the probabilistic model, however, lies mostly on the quality of 
available data. Indeed it is necessary to know the mean value of the significative 
wave height conditional on the wind velocity uy ju = f(U) and the standard 
deviation of H, conditional on the wind velocity U, op ju = f(U). Thus, such 
a stochastic model is strictly recommended when reliable estimates of 4g. ju and 
oH,\u are available. Unfortunately this is not our case. And, in addition to this, in 
the present thesis it is not decisive to have at our disposal real data to implement 
the IFORM because it is not central for the research. 


In similar circumstances, in accordance with [2], the significant wave height H, 
defining the sea state state severity, can be assumed independent on the mean wind 
speed, so that H,,, is found from the marginal distribution of Hs and with the 
same sea state duration as the sea state duration used for the construction of the 
environmental contour. 


Several distribution models for H, and T, are available, [81], [82] [83] [84] [85], 
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[86], [87], [88], among others. Often the Generalized Gamma distribution proposed 
in [12] gives a reliable description. 


Extreme Sea State, simplified definition 


Setting up a probabilistic model to estimate H, and T, (or T,) departing from 
the 50-year return period wind velocity Uso, even by assuming an uncorrelated 
model, would be a mere exercise in this context because we are not aiming at a real 
design. In fact, we are more interested in constructing a valid model that, however 
determined the extreme values for U, Hs and Tp, is able to simulate with a higher 
accuracy the extreme offshore environment and the induced response. 

The simplified approach consists in defining a sea state by assigning a severity 
rank depending on the wind velocity. Each severity class is characterized by a range 
of significant wave height and wave period. 

By virtue of this simplification, the only random variable is the mean wind 
speed U, while H, and T, are deterministically determined from the wind velocity, 
see figure 5.16. Thus, equation (2.4) presented in Chapter 2 becomes 


Prat = | P (My > MylV) PA) (5.1) 


| oe SA 


Figure 5.16: Simplified environmental model: the sea state is defined deterministi- 
cally depending on the mean wind speed U. 


In the present case the mean 50-year return period wind speed at 19.5 m above 
the sea level, obtained by using the profile in figure 5.15, is 31.69 m/s (61.61 Kn) 
which, referring to table 5.5 [10], makes the storm fall in the class of “hurricane-type 
storm”. 

This class is characterized by the following ranges for the significative wave- 
height and mean wave period, respectively 


H, = {21.34 to 35.05m} (5.2) 
T; = {10 to 30s} 


100 
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Sea Description of sea Wind speed Significant Wave period 

states irange wave height range 
Knots ft s 

1 Small wavelets 5-10 0.3-1.4 0.5-5 

2 Large wavelets | 10-14 1.4-3 1-7.5 

3 = Small waves 14-18 3-6 1.4-8.8 

4 Small to moderate waves 18-19 6-7 2.5-10.6 

5 Moderate waves (19-24 7-13 2.8-13.5 

6 Large waves 24-30 13-22 3.8-15.5 

7 Moderate gale 30—40 22-45 4.7-21 

8 Strong gale 40-55 [45-70 6.5-25 

9 Hurricane type storm 55-70 70-115 10-30 


Table 5.5: Definition of sea states according to [10]. 


Since this values are quite general, they could result too conservative and not 
very consistent with those specific for the North Sea, where our turbine is supposed 
to be located. 

Hence, referring again to [10] (table 3.19), specific values for the North Sea 
extreme sea state associated with extreme wind of 37 m/s are 


H, =14m (5.4) 
T, = {15 to 17s} (5.5) 


Moreover, as no consideration about the water depth has been made so far, we 
can refer to the diagram proposed in [10] where for the three exposure levels L-1, 
L-2, L-3 of an offshore platform, the relation between water depth and wave height 
has been found.4 

In our simplified model we just need to collect some realistic values, so if we 
consider our offshore wind turbine like an unmanned platform (L-3), then from 
figure 5.17 it is immediate to see that for a water depth of 20m (approximately 
65.6 ft) the corresponding wave height is about 40ft (12.19m). And for such a 
design condition, the recommended extreme wind velocity and wave period are 
58 Kn (29.84m/s) and 11.6s. 

On the other hand, if we refer to Beaufort scale, see table 5.6, we have that for 
a Wind Force 11 (i.e. wind speed grater than 30.60m/s) and Wind Force 12 (i.e. 


4The three exposure levels are defined as follows: 


e L-1: manned, non-evacuated; 
e L-2: manned, evacuated; 


e L-3: unmanned. 
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Figure 5.17: Water depth dependent wave height for different platform types, [10]. 


wind speed in the range of 29 to 33m/s) the suggested values are, respectively 


H, = 10.25m (5.6) 

Tp = 12.59s (5.7) 
and 

H, = 8.70m (5.8) 

Tp = 11.995 (5.9) 


Note that with respect to table 5.6, the peak spectral period for the JONSWAP 
spectrum is given by Tp = 1.1997; or Tp = 1.28772, where Tı and 7> are some 
characteristic periods. 

Although the above data present some differences, a deeper investigation to 
know which are the most suitable parameters for our case is not necessary. Also 
because in real design context a probabilistic model as described in the preceding 
section should be implemented. 

In literature there are many other proposed tables and diagrams to get sea state 
parameters following the simplified approach. Among others we cite also diagrams 
proposed in [89] and [33]. 


5.4.3 Domain decomposition and breaking waves simulations 
In this paragraph a detailed description of the procedure adopted to simulate 
the coupled wind—waves extreme environment is presented. The procedure has been 


implemented in a computer code which basically executes a certain number of in- 
structions which can be grouped into the following steps: 


1. Environmental Analysis 
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Wave Spectrum Parameter Estimates 


Scale of Wind Speed Open Ocean Areas North Sea Areas 
Beau fort at 19.5 m (Bretschneider) (JONSWAP) 


above sea 


3.50 

3.80 

: at 4.20 
13.5 j. 5. : 4.60 
19. j. j. . 5.10 
24.5 7.2 3.65 2. 5.70 
30.5 È 7. .Z 3. 6.70 
37. 4 7 4.85 7.90 
44. 9. ad j: 8.80 
51. 9.50 
59. 10.00 
>64.0 9.6! .23 10.50 


5 
5 


Table 5.6: Definition of sea states according to Beaufort scale. Table from [14]. 


2. Spectrum choice 

3. Irregular sea generation for a requested space-time domain 
4. Setting the wind turbine location 

5. Free surface elevation and zero—crossing analysis 

6. Solve dispersion relations 

7. Check breaking wave limit 


(a) if NOT: Use standard Morison's equation. 
(b) if YES: Identify all the possible time instants at which waves break and: 
i. Define space-time subdomains 
ii. Invoke the fully nonlinear solver 
iii. Perform fully nonlinear simulations 
iv. Plunging breakers analysis 
v. Get parameters to compute impact loads 
vi. Pass impact loads to FAST 
vii. Perform the fully coupled aero—hydro-elastic analysis 


From the environmental analysis, data about wind and wave conditions should 
be collected. However, in a pure conceptual design phase, it is possible to start 
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from the IEC wind class which provides both the U,ef and the turbulence intensity. 
As already said in the pervious section, in this work a wind turbine in class III is 
assumed. 

The power law wind profile permits to get the mean extreme wind speed at 
19.5m above the still sea level. The sea state correlated to such an extreme wind 
condition can be derived by implementing one of the statistical model described 
above, but as the lack of data justifies, here we use the simplified approach. In fact 
we adopt the sea severity parameters among those recommended by [10]. 

From the environmental analysis we basically get {U,e;, Hs, Tp} and after having 
chosen the most suitable wave spectrum, it is possible to move to phase 3. Here, 
first we consider a space-time domain D (t) = [0, Tsim] X [@min, Tmax]; Where Tsim 
is the total simulation time, £min and rmax are the lower and upper bounds of the 
two-dimensional spatial domain. 

Since the focus is devoted to wind-generated waves and also because we imagine 
that the wind turbine is located in the North Sea, the JONSWAP spectrum is 
adopted. See section 4.1.2 for details. 

By using the spectral formulation it is possible to determine the kinematics, the 
velocity potential and the free surface elevation of every water particle p € D (t). 
The so called inverse approach permits in fact to come back to the temporal domain 
from the frequency domain provided that the phase angle e, lost during the Fourier 
transform, is now assigned in such a way to get a different signal but with the 
identical statistical properties. To this end a uniform normal distributed phase angle 
is assumed. 

The two-dimensional time-depending domain D (t) is reduced by one dimension 
by fixing the wind turbine location xy. For all cases we will have x; = 0 as shown 
in figure 5.14. Moreover, as already mentioned, for all simulations the water depth 
is always d = 20m. 

Subroutine 5 performs a zero-crossing analysis of the free surface elevation at 
Tt 

n(x, t) = 5 an COS (—wnt + En) (5.10) 
n 


where a, is the n-th wave amplitude stemming from the JONWASP spectrum 
an = (/2Snn (7) Awn (5.11) 


The time series 7(x;,t) is analyzed by means of a special subroutine which 
computes all the wave periods and the corresponding wave heights. Firstly all the 
time instants up-crossing the still water level are collected into a vector called 
tup- Then the wave periods can be easily computed as the difference between two 
successive zero up-crossing time instants, while the wave heights are obtained by 
taking the difference between the highest and lowest points of the free surface (crests 
and troughs) within the corresponding period. Hence, for each time series these data 
are stored in two vectors, T and H, respectively. In step 6, Solve dispersion relations, 
as many linear dispersion relations as the dimension of T (or H) are solved in order 
to get the vector of wave numbers k which finally leads to the vector of wave lengths 
L. Step 6 ends with the computation of the steepness vector 


E= (5.12) 


S| i 
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Point 7, Check breaking wave limit, performs a comparison between the steepness 
vector € and the limit steepness vector €. By computing the mean value of H and T 
it is not difficult to verify that the condition of deep water is not always guaranteed, 
thus in establishing the breaking limit the water depth is also taken into account 
by the well known relation [5] 


€, = 0.142 tanh (kd) (5.13) 


0.16 


0.12 


0.1 


(SÌ 
È 0.08 


0.06 


0.04 


0.02 
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0.5 
d/L 
Figure 5.18: Limit breaking steepness ey for different wave lengths and water depth. 


From step 7 there could be two outcomes: either the breaking limit is never 
achieved by any of the steepness in €, or there could be a certain number of waves 
which violate the breaking limit. 

In the first case, which is not the focus here, nothing should be done in particular. 
It is reasonable, in fact, to presume that since no impact will take place, Morison’s 
equation remains valid without any alteration. The benefit stemming from the use 
of the fully nonlinear simulator is anyway not negligible because a more accurate 
wave kinematics can be obtained. 

On the contrary, a different scenario occurs when the breaking limit is violated. 

Referring to figure 5.19, and assuming that the vector € has dimension n, it is 
possible the identify a number nb of times t = {ty1, t03;---;t0,,} at which waves 
will theoretically break. In this way all possible breaking events can be analyzed. 

To visualize the importance of detecting the times at which breaking waves may 
occur we can refer to figure 5.14 and notice that it represents a snapshot taken 
at a generic time when, despite of the heavy storm in action, no wave is breaking 
against the substructure. In contrast, at a generic time #,,, or in particular at ty 
an impact may occur as sketched in figure 5.20. 


max ? 


To perform some preliminary simulations aiming mainly at testing the global 
scheme, we start analyzing only the strongest event among nb possible. In the next 
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Figure 5.19: Example of free surface elevation. All zero up-crossing time instants 
are marked with a black dash and the i-th, with i = 1,... nb, wave period and wave 
height are highlighted red. 


section, when extreme wind and waves will be coupled, all the possible events will 
be taken into account. 

Focusing the attention only on the stronger event, we just need to find the 
maximum value of € and then compute the associated time instant ty... as 


tomax — 1 


tomas = tup(1)+ DI T(R+ Timex )/4 (5.14) 
h=1 


where iba, is the vector index corresponding to the maximum steepness, tup is the 
vector collecting the zero up-crossing time instants, T collects all the wave periods 
in the signal n (x+, t). 

Once the vector t, or the time ty,,,, are known, phase 7b ends and data to be 
passed as input for the fully nonlinear solver can be prepared. 

In fact, in step 7(b)i Define space-time subdomain, neighborhoods of x+ and t3,,.x 
have to be found in order to define a specific space-time subdomain Q (t) C D (t) on 
which Laplace’s equation and the fully nonlinear kinematic and dynamic boundary 
conditions can be numerically solved. 

Note that, whenever not ambiguous, the subscript max will be removed so that 
max Will be simply denoted by tọ. 
The subdomain is defined as follows 


ty 


Q (t) = [ty = ty, ty + ôte] x [x = ÔT, £t + ôx] (5.15) 


where Q (t) has the same meaning of section 4. 
Upper and lower bounds of the sub-domain Q (t) are crucial to set the initial and 
boundary conditions to be passed to the fully nonlinear solver. Caution has to be 
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Extreme irregular sea at t» 


} - Lmin } + Lmax | 


Figure 5.20: Sketch of the wind turbine located at x: = 0 in the 2D spatial domain 
Di = [{min, Tmax] for a time t = tp. 


paid when starting the nonlinear numerical solver by assigning initial and boundary 
conditions derived from the linear theory. To avoid numerical instabilities induced 
by a sudden transition from the linear to nonlinear solution, a spatial ramp function 
is applied in a very short space range. 

The radii ôx, and ĝt, of the neighborhood Q (t) are chosen in order to assure 
a good compromise between the accuracy of the nonlinear wave propagation and 
the computational cost. In particular, it has been found that good values to define 
the sub-domains are twice the mean value of the wavelengths (spatial radius) and 
ranging between 2 and 4 s for the temporal radius 


dx, = 2 x mean (L) (5.16) 
Sty = 2 — 4s (5.17) 


where, as already mentioned, vector L collects all the wave lengths contained in the 
signal n (xz, ty). 

It should be noted, however, that the temporal radius ôt, may significantly be 
increased without the global numerical scheme loses its efficiency. For the sole sake 
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Figure 5.21: The three main models involved in the simulation. Wind: IEC Kaimal 
turbulence model; waves: fully nonlinear Boundary Element Method model; impact: 
analytical model. 


of limiting the computational time, in all the applications presented later on, the 
temporal radius is fixed according to the above range. 


Finally, it is possible to move to point 7(b)ii where the fully nonlinear solver 
is invoked. Before doing this, however, it is necessary to set the initial and the 
boundary conditions which have to be passed to the Boundary Element Method 
solver discussed in section 4.2. 
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Initial conditions 


The spectral initial conditions, namely the initial free surface elevation and the 
velocity potential, are calculated as follows 


N 
>; An, COS ( — Wnti + En) (5.18) 
n=1 
= DIE SIA MOT (5.19) 
> ne d) s dA i i 
where x € Qy = [xy — zi, £e + ôx] and ti = ta — dtp is the initial time of the 


numerical simulation. The subscript t on the domain Q (t) just denotes the space- 
time domain at a given time instant, namely Q; turns out to be a spatial domain. 
In addition to that, as customary, an, kn, Wn, En; denote respectively the wave 
amplitude, circular frequency, wave number and random phase angle associated 
with the n-th wave component. The meaning of these symbols will not be repeated 
in the remaining of the work. 

Moreover, note that in equation (5.19) Wheeler” stretching [68] has been used. 
Accordingly, let q (x) = d/(d+n(x,t;)), then 


n (x) = q(x) n (x, ti) +d (q (£) — 1) (5.20) 


In general, the code enables users to chose the most preferred extension of the 
linear solution up to the actual free surface elevation. Indeed, Wheeler, extrapola- 
tion, constant and no stretching, are all possible switches. 


Boundary conditions 


The transition from the linear to the fully nonlinear solution at both ends T; (t) 
and To (t) of the sub-domain Q (t) is made using a two-sided ramp function Rs, 
which is required to be long not more than 10 times one boundary element length. 
The spatial ramp function is defined as follows 


3 sin (in — 2) + 1 for TE [r ti Lempi] 
Rs (x) = 1 for x € (x; + Lrmp1; Lo — Lymp2) 
3 [sin (gr — 3) +1] for TE [£o — Lrmp2, £o] 


where the shorter notation £i = £t — Ô£t, £o = £t + 62; has been introduced. 

This ramp function is necessary because it should be kept in mind that the we 
are setting up a numerical solver by assigning boundary conditions stemming from 
the linear solution. Thanks to this expedient we can safely assume that the solution 
inside the reduced sub-domain 


Q (t) = [to ui ty, ty + ôte] x [xe "= ÔL + Lingis Tt + ÔL Ga Lrmp2] (5.21) 
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is fully reliable as the effects of the linear boundary conditions are all confined in 
the transition zones Lympi and Lymp2. Note, in particular, that the area of our main 
focus is the very close surrounding of the turbine x+, therefore as it lies always in 
the reduced sub-domain 2 (t), we can totally relay on the accuracy of the numerical 
solution. 

A simple explanatory example on how the ramp function (5.4.3) works is given 
in figure 5.22, where the ramp function Rs is applied to a generic function g = g (x). 


The boundary conditions on T; (t) and To (t) are derived from the linear theory 
as follows 


N 
i cosh (kn (d+ y’ (x 
o (xi, t) = > anWn EA (2) cos (kna; — Wnt + En) (5.22) 


n=1 


N . , 
Ladd) 


sini cos (knto — Wnt + En) (5.23) 


n=1 


where t € Q} = [ty — Oty, te + dt] and y' is given by the general form of equa- 
tion (5.20) 


y (x) =q (x) y (a, ti) +d (q (x) — 1) (5.24) 
with 
y' € [-d, 0] is the computational vertical axis; 
y € [-d,n(x,t)] is the actual vertical coordinate up to the free surface eleva- 
tion; 


q is a dimensionless factor; 


n (x,t) the elevation of the actual free surface elevation up to which the solu- 
tion is sought. 


The subscript x on the domain Q (t) just denotes the space-time domain for a 
given location, so that O, turns out to be a time domain. 

To compute the second-order Lagrangian derivatives, used in the solution of 
the second Boundary Value Problem, refer to section 4.2 for details, also the ac- 
celerations need to be assigned as boundary conditions on T; and Io, therefore we 
have 


N 

di (xi, Y, t) = > AnWn gia (knd) sin (knti — Wnt + En) (5.25) 
= cosh (kn (d + y’)) 

Ù? (£o,y, t) = 2 anw a (kad) sin (kn£o — Wnt + En) (5.26) 


After the initialization, all the kinematic quantities computed by the fully non- 
linear numerical solver at each time step are made compatible with the linear bound- 
ary conditions by using the ramp function over the transition zones in the following 
fashion. 
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(a) Generic function g = g (x) and two-side ramp function Rs. 
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(b) Function g after the application of the ramp function. 


Figure 5.22: Example of application of the space ramp function Rs on a domain 
Q = [-150, 150] with Lrmpi = Lrmp2 = 30. 


Let f” = f” (p,t) be a generic quantity numerically computed by the fully non- 
linear solver, with the pair (p, t) € © (t). For a given time instant t € [ty — ôt», te + ôte] 
and for the points lying on the free surface, i.e. y = yf = n, the function f” (p,t) 
turns out to be depending only on the the x-coordinate and it will be denoted by 
Î" = F” (x) with x € [ex — day, xi + ôx]. Note that f” is used for by, vF, v7 and n 


An integrated nonlinear wind—waves model for offshore wind turbines 111 


itself. 
Moreover, let fe be the respective linear quantity computed analytically. 
Then we have 


f (x) = F° (x) (1 — Rs (£)) + f? (£) Rs (£) Va € [xy dr, + ôx] (5.27) 


Notice that in the present case it is not necessary to have two different transition 
zones so that, from now on, we will always assume Lrmpi = Lrmp2 = Lrmp. See 
figure 5.23. 


Transition zone 


Linear solution 1 1 Fully nonlinear 
TÀ solution 1 
Fully nonlinear 0 AN 0 Linear solution 
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Transition zone 
Linear solution 
| -S Xit Lem l +6 Li-L | 
| - 6x1 + Ox | | 
| - Lmin | + Lmax } 


Figure 5.23: Schematic representation of the transition between the linear and fully 
nonlinear solution. The figure is out of scale. 


To summarize, figures 5.24 and 5.25 show a schematic representation of the the 
main steps discussed above about the the global simulation scheme. 


5.44 Applications 


Before coupling the the fully nonlinear wave kinematic solver with the hydro- 
aero-elastic simulator, some applications are here preliminary presented in order to 
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relations A 
Check break- 
ing wave limit 


(to be continued...) 


Figure 5.24: Diagram of the simulation. Part I 


test the reliability of the numerical scheme summarized in figures 5.24 and 5.25. 
Applications 1 and 2 


In the first example a strong gale is simulated and all the relevant parameters 
are listed in table 5.7. 


Sea severity? “Strong gale” (rank 8) 
Urep, Ho, Tp 37.5m/s, 18m, 12s 
Spectrum Type JONSWAP 
Tsim, Tmin, Tmax 300s, —200 m, 200m 
Turbine Location x 0 


Maximum breaking wave time tẹ 8.208 (but not relevant!) 
dx, 146.59 m 
dt 3s 


Table 5.7: Data relevant to application 1. 


In this case no special time instant tẹ has been requested to start the nonlinear 
numerical solver. In other words a generic time t has been assigned. The purpose is 
only to check the quality of the global scheme for the simulation. 

The second simulation also shows a case of a strong gale and similarly to the 
previous example no special time instant t, is requested to start the nonlinear solver. 
It interesting to notice, indeed, in the subdomain under analysis there could be 
multiple plunging breakers. Data relevant to this application are listed in table 5.8 

Figures 5.26 and 5.27 show the evolution of the fully non linear seas for the 
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Figure 5.25: Diagram of the simulation. Part II. 


Sea severity® “Hurricane type storm” (rank 9) 
Crops Ho, Tp 42.5m/s, 22m, 12s 
Spectrum Type JONSWAP 

Tsim, Zminy Tmax 600s, —200m, 200m 
Turbine Location x, 0 

Maximum breaking wave time ty, 123.58 s (but not relevant!) 
ÔT: 137.24 m 

oth 2s 


Table 5.8: Data relevant to application 2. 


two examples above (tables 5.7 and 5.8, respectively). It is immediate to realize 
the transition from the linear solution to the fully nonlinear one work pretty well. 
Moreover, as the initial conditions have been assigned at a generic time instant, 
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Figure 5.26: Five snapshots of a “Strong gale”. Multiple plunging breakers scenario. 
Red arrows denote the free surface particles velocity and the blue dots the boundary 
element mesh. Input data from table 5.7. 


waves break anywhere into the sub-domains. 
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Application 3 


In the third case we show a simulation having a specified time ty. All parameters 
are listed in table 5.9 


Sea severity’ “Hurricane type storm” (rank 9) 
Uret: Hs, Tp 37.5m/s, 18m, 12s 
Spectrum Type JONSWAP 

Tsim, min; Tmax 600s, —200m, 200m 
Turbine Location x, 0 

Maximum breaking wave time ty, 35.38 s (relevant!) 

ÔT: 128.28 m 

Oty 2s 


Table 5.9: Data relevant to application 3. 


This simulation, see table 5.9 and figure 5.28, is extremely important and shows 
that the model described above works well and meet all the expectations. Such a 
result is actually not trivial because the breaking waves prediction tool elaborated 
above (see scheme in figures 5.24 and 5.25) starts form linear boundary and initial 
conditions to simulate a fully nonlinear event. 

The first subfigure (the upper one) in 5.28 is a snapshot taken at t = ty — ôt», 
namely it shows the free surface configuration when we started the fully nonlinear 
solver. The last subfigure (the lower one) in 5.28 shows the free surface at t = ty+dto, 
that is at the end of the refined simulation or, consistently with the concept of sub- 
domains D, at the end of the “sub-simulation”. So as expected in the middle of the 
temporal sub-domain, i.e. at tẹ (figure in the middle), we see exactly what expected: 
the plunging breaker is “crashing” against the virtual turbine substructure causing 
ad additional load to be carefully considered. From this time on, say up to the 
re-entry, the free surface experience very large curvatures and highest velocity at 
the water jet. This is very well seen in the three lower subplot of figure 5.28. These 
three instants are zoomed-out and reported in figure 5.29. 

A clearer representation of the overturning tongue of water is given in figure 5.30 
where the same of figure 5.29 is shown without velocity vectors and grid markers. 


Application 4 


To test further the stability and reliability of the first part of the numerical tool 
developed in this thesis another application is presented. Simulation parameters are 
listed in table 5.10. The main differences here with respect to results associated 
with data in table 5.7 are a larger wave period which tends to reduce the steepness 
in the average and a longer simulation time. In particular, here we have 7, = 16s 
and Tsim = 600s. Moreover, also in this case, as did in the last two, we require to 
investigate the strongest breaker at x+, that means the time ty is relevant. The free 
surface evolution is presented in figure 5.31. The last three snapshots are zoomed 


116 Coupled wind-fully nonlinear waves model 


Sea severity® “Strong gale” (rank 8) 
Ures, Ho, Tp 37.5m/s, 14m, 16s 
Spectrum Type JONSWAP 
Tsim, Zminy Tmax 1200s, —200 m, 200 m 
Turbine Location x, 0 
Maximum breaking wave time tp 464.15 (relevant!) 
52, 184.34m 

Oty 3s 


Table 5.10: Data relevant to application 4. 


and reported in figure 5.32. 

A clearer representation of the overturning tongue of water is given in figure 5.33 
where the same of figure 5.32 is plotted without velocity vectors and grid markers. 

Also this application confirms the reliability of the global numerical model. From 
the second and third subplots presented in figure 5.31 it is clear that the predicted 
time of breaking, ty = 464.15 s, is satisfactory observed and to compute the maxi- 
mum impact force it is possible to use the most severe combination between wave 
elevation and water velocity among all the configurations around tọ. 


Application 5 


This simulation consist of a particular sea state which is not related to excep- 
tional events. Highly nonlinear events may occur also with moderate wind speed. 
Indeed, by considering the rated wind speed of our wind turbine, that is 11.4 m/s, 
possibly after some reduction because that value is referred to the hub height 90 m 
and not to the conventional value of 19.5 m above the mean sea level, we realize to 
fall in a “Moderate Sea” again according to [10]. This sea severity rank (number 5) 
is characterized by the following ranges 


H, = {2.13 to 3.96 m} (5.28) 
Ty = {2.8 to 13.55} (5.29) 


The wave period range is meant to cover all possible periods over which mea- 
surable energy of the random wave for the particular sea state exists. 

Now, by assuming the median value for the significative wave height and a wave 
period in the lower part of its “energy range” it has been observed that overturning 
breaking waves occur as well. 

Table 5.11 lists all the relevant data for this simulation. Four breaking waves 
are detected in this realization and, as did for the previous cases, we isolate the 
stronger event which is expected to occur at tẹ = 33.99 s. 

We chose a larger sub-domain 500 m long in space and 12.5 min long in time. At 
the location of the wind turbine, i.e. x; = 0, the free surface elevation is represented 
by the time series in figure 5.34. 
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Sea severity? “Moderate waves” (rank 5) 
Ure, Ho, Tp 11.4m/s, 3.05m, 4.758 
Spectrum Type JONSWAP 

Tsim, Tmin, Tmax 900s, —250 m, 250 m 
Turbine Location x 0 

Maximum breaking wave time ty 33.99 s (relevant!) 

dr, 33.99 m 

dt 38 


Table 5.11: Data relevant to application 5. 


The three subplots reported in figure 5.35 show that a very steep and unstable 
wave is forming. The lower one clearly shows the vertical profile of the free surface 
just few instants before overturning. 

Looking at what happens in the second group of snapshots, figure 5.36, the 
classical plunging overturning breaker takes place. However, the wave breaks some 
meters before reaching exactly the turbine location. By comparing this numerical 
experiments with the other presented in this section, it should not lead to wrong 
conclusions. Indeed, in this specific case the resolution is much higher because of the 
dimension of the space sub-domain D;. In this case the numerical domain is 68m 
long (which is 2ôx+), while in the previous cases the domain was in the range of 
200m up to 360m, approximately. Within this scale, few meters of approximation 
are definitely acceptable. 

Anyhow, to judge the quality of this simulation, it is also crucial to stress that we 
are simulating a fully nonlinear phenomenon by departing from boundary and initial 
conditions which have been produced by a the linear spectral theory. Therefore, it 
is reasonable that breaking waves my not occur exactly where and when predicted 
by the linear theory. 

Figure 5.37 shows what already presented in figures 5.35 and 5.36, but by remov- 
ing the free surface particles velocity and the mesh markers a clearer representation 
of the water jet is offered. 


Application 6 


The very last case which deserves some attention is a non—breaking circumstance. 
This additional applications confirms once again the reliability of the numerical tool. 
Table 5.12 collects all key input parameters for this simulation. 

What makes singular this application is that in the “Moderate waves” random 
sea we have generated, none of the harmonics reaches such a high steepness to 
violate the breaking limit. However, we fixed as center for the temporal domain ty 
the time instant at which the maximum steepness is reached. 

In the signal, plotted in figure 5.38, there are 242 wave components with mean 
steepness 0.051, mean wave length 39m, mean wave period 4.94s. 

With the same scheme proposed above to identify tẹ, we proceeded by consid- 
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Sea severity! “Moderate waves” (rank 5) 
Uen HaT, 11.4m/s, 3.05m, 5.755 
Spectrum Type JONSWAP 

Tsim, Zmin, Tmax 1200s, —300 m, 300m 
Turbine Location x; 0 


Maximum breaking wave time tẹ NO breakers are detected! 
dr, 39.62 m 
dt 6s 


Table 5.12: Data relevant to application 6. 


ering the maximum steepness and we found that it is t, = 560.90s. According to 
the general scheme discussed above, we launched the fully nonlinear solver with an 
initial time t; = ty — dty = 557.90s and we ended the simulation at a final time 
given as tp + dt) = 563.908. 

The maximum steepness max (€) = 0.1150, which is below the limit in equa- 
tion (5.13), has been used to start the simulator and as it can be seen form the 
third subplot (from above or below) of figure 5.39 the wave is really very steep 
and asymmetric also with respect to the vertical axis. Moreover, even though there 
is no overturning, a nearly vertical front forms approximately at x+. This kind of 
ambiguous situation, where there is no real breaking phenomenon but the velocity 
of water increases a lot storing a large amount of kinetic energy - which could be 
released on the structure - should be investigated further. In our opinion such kind 
of non-overturning waves should also be counted among the impact events as they 
could also induce impulsive contribution. 

The three subplots in figure 5.39 are shown again in figure 5.40 without the 
boundary mesh and the velocity vectors. 


Some remarks on the numerical solver 


We observe that in general all simulations may stop by obeying two criteria: 
a “natural” one, that is when the simulation time reaches the ty + dt); a “forced” 
one, which applies when the water tongue re-enters into the sea surface. In this 
circumstance, in fact, the irrotational hypothesis on which the entire mathematical 
model is founded vanishes and thus the scheme breaks down. 

The numerical instability of the system at the re-entry is shown in figures 5.41 
and 5.42. 

The zoom-in of the subplots in figure 5.41 are shown in figure 5.42. They clearly 
prove that the simulation is stable up to the contact of the water jet with the free 
surface. 

An additional consideration regards also the utilization of local refinement and 
smoothing subroutines. Notice in fact that a multi-breaking waves scenario may 
occur (e.g. the case shown in figure 5.41). In these circumstances, the refinement 
subroutine should work at the same time on different regions of the boundary. 
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It should be pointed out that all the above results have been obtained by employ- 
ing a number of quadratic boundary elements on the free surface NE; adaptively 
computed as follows 


. Ly 
NE; = ceil | —— = NE} 5.30 
“A (= (L) i ( ) 


where ceil and mean simply compute the larger integer and the mean value of their 
arguments, respectively. NE% is number of quadratic elements per wave length and 
as already said, L collects all the wave lengths contained in n (x = x+, t). The opti- 
mum value for NE% has been found to be included in the range of 30-60. 

Furthermore, it is remarkable the fact that during the simulations above the 
regridding subroutine has been called every five time steps, while the smoothing 
was applied every two time steps. Nevertheless, the code presents an extraordinary 
stability because neither restarting nor local refinements in the regions of the cusps 
were necessary. 

Also surprising is to note that the time step is for nearly all the simulation time 
Tsim kept constant on the value 0.05s. When the water jet is forming however, a 
specific subroutine halves the time step span any time a node of the boundary mesh 
undergoes a displacement larger than 50% of the element’s length. This expedient 
permits to circumvent possible instabilities related to very high velocity gradient. 
See section 4.2 for further details. 
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Figure 5.27: Five snapshots of a “Hurricane type storm”. Multiple plunging breakers 
scenario. Red arrows denote the free surface particles velocity and the blue dots the 
boundary element mesh. Input data from table 5.8. 
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Figure 5.28: Five snapshots of a “Strong gale”. Plunging breaker approximately at 
x: and ty. Red arrows denote the free surface particles velocity and the blue dots the 
boundary element mesh. Input data from table 5.9. 
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Figure 5.29: Plunging breaker: zoom of the lower three subplots of figure 5.28, from 
ty to ty + btp. 


ate 


Figure 5.30: Plunging breaker: zoom of the lower three subplots of figure 5.28, from 
ty to ty + ôte. Free surface profiles alone. 
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Figure 5.31: Five snapshots of a “Strong gale” sea state. Plunging breaker at pre- 
dicted values of x; and ty . Red arrows denote the free surface particles velocity and 


the blue dots the boundary element mesh. Input data from table 5.10. 
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Figure 5.32: Plunging breaker: zoom of the lower three subplots of figure 5.31, from 
ty to ty + dtp. 
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Figure 5.33: Plunging breaker: zoom of the lower three subplots of figure 5.31, from 
ty to ty + ôte. Free surface profiles alone. 
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Figure 5.34: Time series of the free surface elevation at x; = 0 for a “Moderate sea”. 
More input data in table 5.11. 
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Figure 5.35: First three snapshots of a “Moderate waves” sea state. Red arrows 


denote the free surface particles velocity and the blue dots the boundary element 
mesh. Input data from table 5.11. 


126 Coupled wind-fully nonlinear waves model 


BEM Free Surface at t = 38.9315 s 
T T 


7 T 


o 
x 


Figure 5.36: Second three snapshots of a “Moderate waves” sea state. Red arrows 
denote the free surface particles velocity and the blue dots the boundary element 
mesh. Input data from table 5.11. 
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Figure 5.39: Five snapshots of a “Moderate waves” sea state. No breaking waves 
occur and ty is fixed by the maximum steepness. Red arrows denote the free surface 
particles velocity and the blue dots the boundary element mesh. Input data from 
table 5.12. 
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Figure 5.40: The three central snapshots of figure 5.39. Free surface evolution in 
the surrounding of x; for a non—breaking wave case. The entire subdomain is shown 
only in the lower subplot. Input data from table 5.12. 
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Figure 5.41: Breaking down of the numerical scheme due to re-entry of the water 
jet in the sea surface. 


(a) Last time step before breaking down. (b) Re-entry and simulation breaking down. 


Figure 5.42: Details of the overturning spout. Zoom-in of figure 5.41. 
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5.5 Fully nonlinear aero-hydro-elastic coupled model 


In this section the complete fully integrated wind-waves simulation model is 
presented. The main steps and the global numerical framework shown in figures 5.24 
and 5.25 is here generalized and extended in order to include all the possible breaking 
wave events occurring during the extreme sea state. The main implication of this 
feature is that the call to the fully nonlinear solver is made as many times as the 
number nb of plunging breakers detected. 

Moreover, in this section the interface of the numerical wave simulator with 
FAST is developed. 


Spectra choice Hazard 


cons» | CSO) Analysis 


M, Kaimal, etc.) 


Simulation domain 
Irregular sea Generate 
Dr > 5 <— 
generation random phase € 
[E min; Emax] X [0, Tsim] | 
Í Solve disp. relations Zero—crossing Fix WT location 
2 = gk;tanhk;d | is i u 
wy = gk; tanh kj analysis in x Li € [min, Tmax] 


Check breaking <> compared with 
wave limit (H/L) 


breaking 


(U, Tp, Hs) 


(to be continued...) 
Figure 5.43: First part of the complete aero-hydro-elastic simulation. 


The diagram shown in figure 5.43 does not present any significant difference 
compared with the scheme already seen in figure 5.24. It is just repeated for the 
sake of completeness. 

On the contrary, a substantial difference appears in the second part of the com- 
plete simulation scheme depicted in figure 5.44. Here in fact all the instants at which 
a breaking wave is expected are collected in the vector tẹ and nb sub-domains have 
to be defined. In fact, while in the previous section only the strongest event was 
considered, here we need to account for all the possible impacts which may occur 
during a storm. Therefore, nb possible overturning breakers are simulated by call 
ing nb times the fully nonlinear solver. Note that the procedure shown earlier to 
compute the time instant at which braking occur needs to be slightly adjusted in 
order to get the vector ty 


BrLm(i)-1 
to; =tup(1)+ XO T(h) +7 (Brim (i) /4 (5.31) 
h=1 


where BrLm is the vector selecting only elements of € which breaks, tup is the vector 
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Figure 5.44: Second part of the complete aero-hydro-elastic simulation. 
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Figure 5.45: Definition of subdomains and initial and boundary conditions assign- 
ment form JONSWAP spectrum. 
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Figure 5.46: Last part of the global simulation scheme: interface with FAST. 


The definition of sub-domains is well depicted in figure 5.45 where it is clear how 
the BEM-based code is activated to solve the fully nonlinear Laplace’s equation on 
the domain Q; (relative to the i-th braking event) by assigning initial and boundary 
conditions derived by the linear spectral approach. 

The last group of subroutines of the global simulation code execute the in- 
structions sketched in figure 5.46, where the kinematic properties of each plunging 
breaker are derived. Impact forces are computed and subsequently passed to FAST 
which performs the time marching analysis. 

The instructions executed by the software are also summarized in the algo- 
rithm 3 which gives a simplified representation of the whole procedure. Note that 
the environmental analysis is carried out externally from algorithm 3. 

In the above algorithm WaveDT denotes the time step for the irregular sea 
simulation, the argument input; passed to the subroutine CallBEMSolver collects 
all the data related to the ith breaking wave necessary to run the fully nonlinear 
code. See section 4.2 for details. 

A global view of all the main steps of the code is presented in the scheme on 
figure 5.47. 

Before presenting the applications of the software, we just stress that fact that in 
the light of the considerations and values proposed in the previous section about the 
simplified extreme sea states definition, given a wind turbine located in a generic re- 
gion of North Sea with water depth of 20 m designed to withstand an 50-year return 
period extreme mean wind velocity of 37.5 m/s (wind class TIT, [36]), it seems to be 
reasonable to perform most of the simulations by assuming variations of the input 
parameters in the neighborhoods of the “most probable” sea state characterized by 
the following values of significative wave height and peak period 


e H,=12m 
eT, = 12s 


Now it is possible to start simulating irregular seas coupled with realistic tur- 
bulent wind. 
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Figure 5.47: Framework of the whole simulation scheme. 
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Algorithm 3: Global simulation scheme. 
input : Environmental variobles: Uref, Hs, Tp 
input : Total simulation time: Tsim 
output: Impact force: fr 
Initialization: fr = zeros (0 : WaveDT : Tsim, 3); 


Spectrum choice; 

Irregular sea generation; 

Setting the wind turbine location; 

Zero-crossing analysis; 

Solve dispersion relations; 

Check breaking waves limit; 

if breaking waves occur then 

Identify all possible breakers: get nb; 

for i = 1 to nb do 
Define space-time subdomains: fix dtp, and 624; 
[vi, M,| = CallBEMSolver (input;); 
[fr,] = ComputelmpactForce (vi, m, ); 


else 
Morison’s eq. unaltered; 


Get final impact force fp = 37”, fi; 


5.5.1 Fully coupled Simulation #01 


Before approaching the final coupled simulations, a preliminary model is set up 
with no wind blowing, while an extreme irregular sea is generated. We should point 
out that this scenario is unrealistic because, as stated many times throughout this 
work, the focus is on wind-generated waves and not swells. However this preliminary 
simulation allows to test whether the software, until now only tested with regular 
waves, is stable and reliable. 

In table 5.13 all the input data characterizing this simulation, form now on 
denoted by Simulation #01, are provided. 

Note that WaveDT, DT, dt in the above table refer to time step for the linear 
irregular sea generation, FAST time step, our BEM simulator, respectively. More- 
over, values dx; and ôt, will always be constant for each breaking event occurring 
in one simulation. 

First of all we start by showing the time series of the free surface elevation at 
the turbine location. The first time the breaking wave limit is achieved is at 68.12 s 
when, as confirmed in figure 5.48, a clear peak of about 13 m is registered. 

Since nb = 2, the fully nonlinear wave kinematic solver is called two times. The 
most significative snapshots of the free surface evolution for each call are shown in 
figures 5.49 and 5.51, respectively. 

Data passed to the subroutine computing the impact force are extracted at the 
time shown in figure 5.50. 

It is interesting to observe that the predicted breaking wave time matches very 
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Uref, Iref 0 m/s, 0 
U19.5, Hs, Tp 0m/s, 12m, 12s 
Spectrum type, y JONSWAP, 3.3 


Water depth 


Intermediate Water 


Tsim; Tmin, Tmax 


600s, —300 m, 300m 


WaveDT, DT, dt!! 


0.05s, 0.01s, 0.05s 


Turbine location x; 0 
Predicted number of breaking wave events 2 
Number of breaking wave events occurred nb 2 
Predicted breaking waves time vector tp s 68.118, 105.376 
BrLm 8, 11 

Oat 252.7681 m 
dt 3s 

NE; 80 


Table 5.13: Input data for Simulation #01. 
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Figure 5.48: Time series of the free surface elevation at x; = 0 for Simulation #01. 


Input data in table 5.13. 


well the instant at which the impact occurs. However, it will not always be like this. 


It must be pointed out that the predicted time at which the breaking wave 
would occur is just an indication to define the initial time to start the BEM solver. 
In fact, as it happens for the second breaker, the predicted time was 105.38s but 
the simulation reveals that the wave crest starts overturning (say when the front 
forms) later at t = 107.70s. From figure 5.53, it is clear that the overturning is 
taking place approximately 20 m further the turbine location. 

This example reveals how inaccurate the linear model can be in predicting break- 
ing waves. 

The total impact force vector read in FAST is shown in figure 5.55. This plot 
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Figure 5.49: Simulation #01, first 
nonlinear free surface evolution. 
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Figure 5.50: Simulation #01, first breaking wave event. Time at which the slam is 
supposed to happen. 


is rather useful in understanding how the impulsive contributions are distributed 
during the whole simulation. 


The structural response is presented in figures 5.56 and 5.57. The two slamming 
loads which occur during 10min simulation induce on the structure a remarkable 
additional state of internal stress. To have a proper idea of what really happens when 
these two impulsive actions act on the monopile, we first observe the time series of 
the tower top fore-aft displacements, figure 5.56. From a maximum value (either 
positive or negative) of 5cm when the first impact takes place, the maximum peak 
displacement leaps to approximately 21cm, then when the vibrations start decay 
suddenly the second hit arrives and, even though it has a lower intensity because of 
a smaller m and reduced impact velocity, a clear amplification is registered around 
105s which brings the peak displacement up to 24cm. 


Note that here ultimate strength conditions are investigated, therefore the peak 
values are really crucial in assessing the structural safety. 


Also looking at the shear force and bending moment at the foot of the monopile, 
reported in figure 5.57, the two peaks associated with the respective plunging break- 
ers are very clear. In this case, contrarily to what happens for the tower top fore-aft 
displacement, the second shock does not cause an amplification of the first one. Both 
Fx and My, indeed, seem to have the time to dissipate all the momentum induced 
by the first impact in about 35s which is approximately the time interval between 
the two slams. 


The remaining internal forces and moments are also slightly affected by the 
impacts, but of course, due to the main direction of the colliding waves, they register 
a negligible consequences and for this reason are not shown in this first simulation. 


From this preliminary case analyzed interesting suggestions to set up the next 
case study are derived: first the simulation time should be increased, second it is 
fundamental now to activate the turbulent wind simulator in order to have a first 
more realistic extreme environment. 
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Figure 5.51: Simulation #01, second breaking wave event. Four snapshots of fully 
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Figure 5.52: Simulation #01, impact forces computed according to section 4.5. 
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Figure 5.53: Simulation #01, second breaking event. Time at which the slam is 
supposed to happen. 
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Figure 5.54: Simulation #01, second breaking event. Time at which the slam is 
supposed to happen. Detailed view of the impact front of the jet forming shown in 
figure 5.53. 
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Figure 5.55: Simulation #01, time history of the impact forces throughout the total 
simulation time. On this scale the two impacts look like just two pins with intensity 
in agreement with figure 5.52. 
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Figure 5.56: Simulation #01, Tower top fore—aft displacements time series. 


5.5.2 Fully coupled Simulation #02 


Table 5.14 lists all the input data characterizing this second simulation, form 
now on referred to as Simulation #02. This case is much more realistic than the 
previous one because here we really have for the first time the coupled wind and 
wind-waves extreme actions together. 

The turbulent wind is generated by means of TurbSim, already introduced in 
section 5.4 of this chapter. The three 50-year recurrence period turbulent wind 
velocity components are plotted in figure 5.58. 

Such a severe wind condition gives rise to an extreme irregular sea modeled by a 
JONSWAP spectrum characterized by a significative wave height Hs = 12m and a 
peak period T, = 12s as reported in table 5.14. This sea has a free surface elevation 
at the turbine location shown in figure 5.59. 

The software developed in this thesis at this point performs a zero-crossing 
analysis from which it results that in 10 min simulation and upon the severity en- 
vironmental conditions assumed, one breaking wave event may occur, i.e nb = 1. 

Thus the fully nonlinear wave kinematic solver is called only once. The most, 
significative snapshots of the free surface evolution for the call are shown in fig- 
ure 5.60. 

Figure 5.60 shows four frames of the simulation. The upper subplot represents 
the initial time t;, see equation (5.15), at which the BEM solver is started. The two 
intermediate sub-figures show the augmentation of the steepness until the “cusp” 
starts overturning, while the lower one describes the re-entry of the water jet which 
is completely curled. The latter is also the last time step the software in able to 
integrate the boundary conditions, in fact after this, the multi-connected domain 
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Figure 5.57: Simulation #01, tower base shear force Fyt and overturning moment 
My. 


occurs, the irrotational flow hypothesis is no longer valid and the numerical scheme 
breaks down as shown in the example of figure 5.41. 


Data passed to the subroutine that computes the impact force are extracted 
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Uref, Iref 37.5 m/s, 0.12 
Ui9.5, Hs, Tp 31.69 m/s, 12m, 12s 
Spectrum type, y JONSWAP, 3.3 


Water depth 


Intermediate Water 


Tsim; Tmin, Tmax 


600s, —300 m, 300m 


WaveDT, DT, dt!? 


0.05s, 0.01s, 0.05s 


Turbine location x; 0 
Predicted number of breaking wave events 1 
Number of breaking wave events occurred nb 1 
Predicted breaking waves time vector tp s 516.10 
BrLm, Steepness e 48, 0.1473 
ÔT: 263.304 m 
dt 3s 
NE; 100 


Table 5.14: Input data for Simulation #02. 
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Figure 5.58: Time series of three turbulent wind speed components according to 
the Extreme Wind speed Model (EWM) of IEC61400-1 3rd ed. Time histories used 


in Simulation #02. 
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Figure 5.59: Time series of the free surface elevation at x; = 0 for Simulation #02. 
Input data in table 5.14. 


when the worst condition among maximum elevation and maximum impact velocity 
is reached. This condition is shown in figure 5.61. 

Also in this circumstance the predicted breaking wave time matches extremely 
well the instant at which maximum wave elevation and impact velocity are ex- 
tracted. The predicted value is 516.10s, while the time at which the crest starts 
overturning and may likely exert the strongest impact force is t = 516.85 s. The two 
value are surprisingly close each other and this confirms the validity of the approach 
hereby implemented. A zoom-in of the free surface elevation with its kinematic fea- 
tures is offered in figure 5.62. 

Note also that the spatial agreement between the expected location a, = 0 
(where the turbine is supposed) and the simulated impact. 

The elevation m and the velocity v characterizing this impact are 12.44 m/s and 
7.73m, respectively; the impact duration is 0.0978s and this means that approx- 
imately 10 time steps of FAST integration scheme account for such an impulsive 
contribution. 

The total impact force vector read-in by FAST is plotted in figure 5.63(b). 
This plot permits to better understand the following results about the structural 
response. Figures 5.64 shows the tower top fore-aft displacement time series and, 
contrarily to what happened in Simulation #01, in this case impulsive load does 
not bring the tower top displacement too mach beyond its normal maxima. Where 
with “normal maxima” is intended as the peaks in the displacement time series 
before tp. 

The reason why this happens lies in the fact that Simulation #01 was performed 
without wind and, of course, this made the structural displacement only depend 
on the the hydrodynamic loads where the impulsive action was generated. On the 
contrary, in Simulation #2 a turbulent extreme wind (with mean value of 37.5 m/s) 
is blowing in the same direction of the slamming wave. Therefore, the tower is 
already vibrating with an average maximum displacement of 20cm (it was 5cm 
for Simulation #01) and when the impact occurs it can even be out of phase. We 
want to point out that the frequency of the impulsive loads are very high, therefore 
they are far beyond the first system fore—aft natural frequency. For this reason, the 
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Figure 5.60: Simulation #02, first breaking wave event. Four snapshots of fully 
nonlinear free surface evolution. 
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Figure 5.61: Simulation #02, maximum wave elevation n, and impact velocity v 
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Figure 5.62: Simulation #02, second breaking wave event. Time at which the slam 
is supposed to happen. Detailed view of the impact front of the jet forming shown in 
figure 5.61. 


consequences of the impacts are almost negligible on the global dynamics of the 
system. 

The slamming wave which may occur during 10 min simulation induces on the 
structure a remarkable additional state of internal stress. And likewise to what 
pointed out for Simulation #01, peaks of internal forces and moments are abso- 
lutely relevant for our scopes. 

It is impressive but not surprising to see the instantaneous augmentation of both 
the shear tower base force F,, and the bending tower base moment My; plotted 
in figure 5.65. The force undergoes a leap from a peak of 4 x 10° kN to nearly 
10 x 10 kN (positive sign), while the bending moment has extreme peaks of about 
0.5 x 10° kNm and the impact brings the value up to nearly 1.5 x 10° kNm. 

Another interesting structural output which deserves attention is the internal 
bending moment at the still water level. In this cross-section the structure registers 
the additional couple due to the impact given by the impact force My = Frdr, 
where Fr = frAm is the impact force described in section 4.5 and dr denotes the 
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(b) Time history of the impact force throughout the total simulation time. On this scale the impact 
looks just like a spike. 


Figure 5.63: Simulation #02, impact force time histories in the two different time 
scales. 


moment arm given as 7 (1 —1/2A), where A is the curling factor [8]. 


Also for the fore-aft bending moment at the tower cross-section taken at the 
mean sea level, plotted in figure 5.66, as expected, a peak is registered at t = 517.5 s 
due to the impact. Note also that since in this case the moment arm is twenty meters 
(the water depth) shorter with respect to the moment acting at the monopile foot 


form plot 5.65, M,; results one order higher: 1.5 x 10° kN versus 3.4 x 10° kNm at 
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Figure 5.64: Simulation #02, Tower top fore—aft displacement time series. 


the SWL!3 

Since for this simulation lateral and vertical wind components are also simulated, 
see figure 5.58, it would be interesting also to see the structural response in these 
other directions, to check whether the impact causes some effects. 

As expected, forces, bending moments and torsion at the tower base are not com- 
pletely negligible due to the lateral and vertical wind loads, however from figure 5.67 
it is possible to say that no relevant effects are induced by the impact. 


5.5.3 Fully coupled Simulation #03 


The simulation presented here is particularly important for two main reasons: 
first it is longer than those launched up to now, and this is fundamental to have 
much more reliable data when evaluating the distribution of peaks due to impacts. 
The duration is also relevant because it permits to test the numerical stability of 
the code. This application shows that the code developed is reliable also for 40 min 
simulations. 

The second very important reason for which Simulation #03 is crucial is that 
it is based on “real data” input. As discussed in section 5.4, so far we have assumed 
the most important data to estimate the sea severity by using a simplified approach. 
We just relied on data and tables provided in literature and said that for 31.695 m/s 
wind speed at 19.5m above the sea level, the related sea state severity should be 
defined in a neighborhood of H, = 12m and 7, = 12s. This is true but in this 


13 This numbers are totally indicative and have the only scope to make rough check whether the 
results match the expectations. 
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Figure 5.65: Simulation #02, tower base shear force F,t and overturning moment 
Myt. 
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Figure 5.66: Simulation #02, tower top fore-aft bending moment at the tower 
cross-section taken approximately at the mean sea level. 


case we wanted to use input parameters which have been extrapolated by real data 
referred to the North Sea. In [31], indeed, starting from data recorded over 3 years in 
the area of 2° longitude east and between 53° and 54latitude north, the significant 
wave height and mean zero-crossing period with return period 50 years have been 
extrapolated. So in this case we are using data that we are sure they have the same 
return period of the wind velocity assumed. Anyhow, it should be kept in mind that 
the distributions for H, and 7, have been assumed uncorrelated to the mean wind 
speed and this could lead to some approximation of the values. All the parameters 
characterizing Simulation #03 are listed in table 5.15. 

The turbulent wind is generated as in the case of Simulation #02 and the 
three 50-year recurrence period turbulent wind velocity components are plotted in 
figure 5.68. 

The 50-year return period random sea which may occur in the North Sea during 
a storm having the turbulent wind plotted in figure 5.68 is depicted in figure 5.69. 

The zero-crossing analysis of the signal plotted in figure 5.69 detects that there 
could be 3 wave components with such a high steepness to break. Hence, the fully 
nonlinear solver is invoked three times (nb = 3). For each of these sub-simulations, 
three significant configurations of the free surface evolution are shown in the follow- 
ing. In all three cases, the upper subplot represents the initial time when the BEM 
solver is called. The three breaking wave events are shown in figures 5.70 to 5.72 

The last expected impact would occur at t = 470.03s. And according with the 
sub-domain defined above, the fully nonlinear simulator is started with an initial 
time t; = 467.03 s. In this case the steepness is e = € ((BrLm (3)) = 0.096 which, for 
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Time Series Plots of Tower base moments Mxt and Mzt 


x 105 
1 T T T T T T 
E 
zZ 
= 
È 
v 
Q 
H 
z 
E 
È 1 L 1 1 1 fi 
0 100 200 300 400 500 600 700 
Time (sec) 
5000 
E 
z 
= 
v 0 
N 
= \ 
0 
a 
4 
= 
B 
-5000 
0 100 200 300 400 500 600 700 


Time (sec) 


Figure 5.67: Simulation #02, tower base shear force F,, and overturning moment 
Myt. 
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Uref, Iref 37.5m/s, 0.12 
U19.5; Hs, Tp 31.69 m/s, 9.80m, 13.63s 
Spectrum type, 7 JONSWAP, 3.3 
Water depth Intermediate Water 
Tieni 0 mas 2400s, —300 m, 300m 
WaveDT, DT, dt!4 0.05s, 0.01s, 0.05s 
Turbine location x; 0 

Predicted number of breaking wave events 3 

Number of breaking wave events occurred nb 2 

Predicted breaking waves time vector tp s 22.28, 36.05, 470.02 
BrLm, Steepness e 2,3, 41 

ÔT, 286.27 m 

dt 38 

NE; 80 


Table 5.15: Data relevant to Simulation #03. 


a case of intermediate water, is not too larger than the breaking limit. Moreover, to 
increase the numerical stability of the code the smoothing and regridding of the free 
surface have been used frequently: smoothing at each time step, regridding every 
three steps. This causes an unwanted dissipation of energy, which may delay the 
formation of the plunging breaker. 

During the third simulation the expected overturning wave does not occur. See 
figure 5.72. However, a very steep wave forms and its effect cannot be neglected 
when it approaches the monopile. Although much weaker, an impact can also be 
induced. This last case shows that the prediction of breaking waves based on the 
linear theory may result inaccurate. 

The three impact forces stemming from the above simulations are shown in 
figure 5.73. 

The distribution of the impact forces over the whole simulation time is shown in 
figure 5.74. This plot is rather useful to understand the following results about he 
structural response. From this plot two remarks rise: the first two impacts are very 
close each other and this could induce some amplification in the global dynamics 
of the structure. Second, the third steep wave, which does not really plunge, is also 
considered an impulsive action. However, its associated impact force is much less 
strong due to the fact that without plunging, the crest does not reach those very 
high velocities of the other cases. 

Figure 5.75 shows the tower top fore-aft displacement time series. It seems 
that the two impacts produce an increase in the displacement with a certain delay. 
Consider that the structure is undergoing a strong turbulent wind and its natural 
frequency is very far away from the impulsive action. 

Concerning the state of stress, the two slamming waves that occur during the 
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Figure 5.68: Time series of three turbulent wind speed components according to 
the Extreme Wind speed Model (EWM) of IEC61400-1 3rd ed. Time histories used 


in Simulation #03. 
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Figure 5.69: Time series of the free surface elevation at x; = 0 for Simulation #03. 
Input data in table 5.14. 


40 min long simulation induce on the structure a remarkable additional state of 
internal stress. 

From output shown in plots 5.75 and 5.76 two considerations can be directly 
done. First of all also in this case the instantaneous augmentations of both the 
shear tower base force Fy; and the bending tower base moment M, are phenome- 
nal. They bring the average peak value from approximately 2 x 10° kNm to nearly 
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Figure 5.70: Simulation #03, first breaking wave event. Three snapshots of fully 
nonlinear free surface evolution. 


8 x 10° kNm. The second thing that can be observed is the effect of the third non- 
breaking wave. Since it has included into the simulation as a breaking wave, the 
force it exerts is almost in the range of the standard maxima due to plain Morison’s 
equation. Therefore, this instructive simulation tells us that without the formation 
of the jet the effect on the structure is totally different. The last remark we can add 
is that the peaks tend to decay very rapidly, indeed the negative peak is already 
less intense. This of course much depends on the global system dynamics. 


Another interesting structural output which deserves attention is the internal 
bending moment at the still water level. In this cross-section the structure registers 
the additional couple due to the impact given by the impact force My = Frdr, 
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Figure 5.71: Simulation #03, second breaking wave event. Three snapshots of fully 
nonlinear free surface evolution. 


where Fr = frAm is the impact force described in section 4.5 and d; denotes the 
moment arm. 


5.5.4 Fully 


coupled Simulation #04 


This simulation is characterized by the data set reported in table 5.16. The wind 
conditions are unchanged with respect to the previous case, while sea state severity 
parameters have been slightly modified. The total simulation time is 15 min. This 
simulation differs from the others mainly in the number of impacts. 


The three components of the extreme turbulent wind are plotted in figure 5.78. 
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Figure 5.72: Simulation #03, third (expected) breaking wave event. Three snap- 
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shots of fully nonlinear free surface evolution. 


The zero-crossing analysis of the free surface at the wind turbine location detects 
four breaking waves, thus the fully nonlinear plunging breaker simulator is called 


four times, i.e. nb = 4. 


For each of these sub-simulations three significant configurations of the free 
surface evolution are shown in figures 5.79 to 5.82. In all four cases, the upper 


subplot represents the initial time when the BEM solver is called. 


From figures 5.79 to 5.82 it can be seen that the first breaking wave event is not a 
real overturning phenomenon. Or at least the simulation duration [ty, — ôte, to, + ôte] 
was not enough to capture the complete plunging. There are many reasons for which 
sometimes the predicted breaking event does not actually happens. For sure it must 
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Figure 5.73: Simulation #03, impact forces computed according to section 4.5. 
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Figure 5.74: Simulation #03, impact forces time history over the all simulation 


time Tsim. 
Usa Iref 37.5m/s, 0.12 
Uiss, He, Tp 31.69m/s, 11.50m, 10.60s 
Spectrum type, 7 JONSWAP, 3.3 


Water depth 


Intermediate Water 


Tsim; Tmin, Tmax 


900s, —300 m, 300 m 


WaveDT, DT, dt!5 


0.05s, 0.01s, 0.05s 


Turbine location x; 0 

Predicted number of breaking waves 4 

Number of breaking waves occurred nb 3 

Predicted breaking waves times fp s 72.15, 114.56, 482.49, 551.06 
BrLm 7, 11, 53, 60 
Steepness e (at breaking) 0.134, 0.133, 0.146, 0.144 
dr: 205.02 m 

dt 2s 

NE; 80 


Table 5.16: Data relevant to Simulation #04. 
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Figure 5.75: Simulation #03, tower top fore—aft displacement time series. 


be remarked that the breaking limit criterion is a theoretical value which in certain 
situations (e.g. linear wave theory, as in our case) can be less reliable, secondly it 
could happen that the very small amount of energy dissipated by the smoothing 
subroutine, which is generally not affecting the solution, in this sensitive cases can 
influence. 

However, the crest shown in figure 5.79 is steep enough to think that when 
hitting a pile a certain impulsive effect is also produced. Therefore it has also been 
included in to the total number of impact forces passed to FAST. 

The four impact forces stemming from the above simulations are shown in fig- 
ure 5.83. 

The distribution of the all impact forces over the whole simulation time is shown 
in figure 5.84. As for the previous case, this plot is rather useful for qualitative 
judgment of the extreme structural response. From this plot it is evident that the 
impact associated with the third breaking wave, represented in figure 5.81, is the 
strongest. Even though its elevation is not higher than the other cases, the water 
particles at the cusp have an average velocity of 15.05 m/s which makes the crash 
really strong. 

Figure 5.85 shows the tower top fore—aft displacement time series. It seems that 
the first two impacts induce an increase of the tower top displacements with peaks 
of about 25cm. The third impact does not produce an instantaneous effect on the 
displacement, anyway amplifications are registered during the third a fourth slams. 

On the contrary, the tower base forces and moments reflect instantaneously 
the effects of the impacts. Peaks in figure 5.86, in fact, strictly follow the impacts 
distribution over the all simulated time as shown in figure 5.84. 
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Figure 5.76: Simulation #03, tower base shear force Fyt and overturning moment 
Myt. 
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Figure 5.77: Simulation #03, tower top fore-aft bending moment at the tower 
cross-section taken approximately at the mean sea level. 
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Figure 5.78: Time series of three turbulent wind speed components according to 
the Extreme Wind speed Model (EWM) of IEC61400-1 3rd ed. Time histories used 


in Simulation #04. 
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(c) Time at which data to compute the impact force are taken. 


Figure 5.79: Simulation #04, first breaking event. Three snapshots of fully nonlin- 
ear free surface evolution. 
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(b) Generic intermediate time step. 
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(c) Time at which data to compute the impact force are taken. 


Figure 5.80: Simulation #04, second breaking event. Three snapshots of fully non- 
linear free surface evolution. 
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Figure 5.81: Simulation #04, third breaking event. Three snapshots of fully non- 
linear free surface evolution. 
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(c) Time at which data to compute the impact force are taken. 
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Figure 5.82: Simulation #04, forth breaking event. Three snapshots of fully non- 


linear free surface evolution. 
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Figure 5.83: Simulation #04, the four impact forces computed according to sec- 


tion 4.5. 
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Figure 5.84: Simulation #04, impact forces time history over the all simulation 
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Figure 5.85: Simulation #04, tower top fore—aft displacement time series. 
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Figure 5.86: Simulation #04, tower base shear force Fyt and overturning moment 
Myt. 


Chapter6 
Achievements and final remarks 


This concluding chapter summarizes the principal achievements of this thesis. More- 
over, the important implications of the proposed numerical tool on the structural risk as- 
sessment phase are underlined. 


6.1 Implications of the proposed model on the structural safety and Risk Assess- 
ment 


At this point of the work some crucial and spontaneous questions rise up: what 
are the practical consequences of such a model? Why should it be used instead of the 
traditional impact models? Would not it be easier to use directly the semi-empirical 
deterministic impact wave models available in literature? 

Let us start from the last issue. Modern design approaches aims at employing 
tools (especially numerical tools in this context), which are more and more rep- 
resentative of the reality (reduction of model uncertainties) and at the same time 
guarantee an acceptable cost in terms of computational time. In fact, to minimize 
the model uncertainties, one could model in a very refined scale the fluid—structure 
interaction. Unfortunately, this approach does not meet the second fundamental 
need when designing offshore wind turbine, that is the limited computational time 
available. 

The new numerical model developed in this thesis precisely represents a com- 
promise for the two major needs above mentioned. It is able, indeed, to simulate 
offshore monopile-supported wind turbines exposed to wave impacts with very high 
accuracy without penalizing the computational effort normally required. This is pos- 
sible thanks to two main features of the computational strategy developed: i) the 
domain decomposition linear-nonlinear; ii) the implementation of the Boundary 
Element Method-based solver which, among others, has the brilliant advantage to 
reduce the problem under consideration by one dimension. 

Usually, impact forces are computed with some semi-empirical formula derived 
mainly from experiments. Just to give an example, in [12] and [10] it is recommend 
to compute the impact pressure p as follows 


p = priv? (6.1) 


where «kı = 5.98 for impact due to waves breaking in proximity of the vertical 
cylinder (this is always our case), while kı = 2.74 when the pile is hit by broken 
waves. p is the water density and v is the impact velocity given as follows 


T 
v = 0.48 x 1.0922" (6.2) 
2T 
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Figure 6.1: Impact force through its duration taken from Simulation #01. 


If we recall the impact in Simulation #02, it is possible to compare the the sim- 
ulated impact velocity with the predicted one by the above semi-empirical model. 
The simulated impact velocity was v = 12.44m/s, while the wave period associ- 
ated with that wave component was T = 10s. With these data the above model 
would lead to an impact velocity of 8.14m/s, while our simulation shows that the 
strongest impact occur when the velocity reaches 12.44m/s. Figure 6.1 is recalled 
here to remind some key values concerning this impact. 

The difference between the two impact velocities gives an idea on how large is 
the range of variation and approximation of semi-empirical models. Therefore, the 
first thing to say about the model developed in the thesis is that it is able to capture 
the wave impacts much more accurately and thus forces transferred to the structure 
are much more reliable. This, by the way, is what we announced in Chapter 1 when 
the reduction of model uncertainties was mentioned. In fact, by using our model 
the extreme response of offshore wind turbines can be more realistically captured 
and this induces remarkable consequences on the ultimate loads assessment. 

The implications mentioned above bring us back to the issue opened in Chapter 2 
when the crucial phase of the Risk Assessment was discussed. As we said, here the 
scope is the quantification of structural risk. To this aim, by recalling the definition 


Risk = Probability of failure x Losses [Losses unit /time] 


it appears that the prediction of the state of failure is deeply conditioned on the 
peaks distribution of each simulation, therefore the utilization of the proposed nu- 
merical approach becomes of primary importance at this stage of the general risk 
management framework. 

In this paragraph the important consequences of the proposed model on the 
vulnerability term of equation (5.1), written as 


Vulnerability = P (My: > My lU = Uso) (6.3) 
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are highlighted by comparing the structural response obtained with the proposed 
wave impact model with the usual non-impact linear wave approach. 

Notice from equation (6.3) that the approach followed in this thesis to investigate 
the extreme response of offshore wind turbines is the one referred to as “external 
conditions—based” approach, see Chapter 2. This means that the extreme response 
is not extrapolated from simulations at wind speeds in the range of the cut—in and 
cut-out, but directly reproduced by assigning extreme loads. 

After running a certain number of simulations, whose input (loads) differ each 
other only by the seeds used to generate the extreme wind and waves, statistical 
properties of each realization should be analyzed. Together with the randomness 
of the sea and the turbulent wind, to make even wider the difference between two 
different arbitrarily selected simulations is the randomness of the impulsive contri- 
butions due to breaking waves. The latter makes necessary to have at disposal a 
minimum number of runs. Due to limited time, here it is not possible to run a large 
number of simulations at least 10 min long each. Hence, just to describe the proce- 
dure it will be assumed that one selected time series for the case considering the 
impacts forces and one realization without impact forces are selected and assumed 
to be representative of all the others. Moreover, for all the remaining of the work, 
the structural response taken as reference is the overturning moment at the tower 
base My. 

Let us consider the case of Mys associated with Simulation #04 which is recalled 
in figure 6.2. Subplot 6.2(a) represents the response when no impact is considered. 
Let us say that this is the output that one would obtain without the numerical tool 
developed in this thesis. On the contrary, subplot 6.2(b) shows for the same sea 
severity, established by the same extreme turbulent wind, the response including 


the wave impacts. 
The statistics of the series in subplot 6.2(a) is given as follows 


Minimum Mean Maximum StdDev Skewness Range 
-8.291e+004 5.198e+003 8.538e+004 1.893e+004 -2.272e-003 1.683e+005 


while the statistics of the series in 6.2(b) is the following 


Minimum Mean Maximum StdDev Skewness Range 
-9.632e+004 5.256e+003 1.601e+005 1.979e+004 3.494e-001 2.564e+005 


A much clearer comparison is possible by the superposition of the two responses 
as plotted in figure 6.3. To extract local maxima form the series a peak—over— 
threshold is used [90], [91]. The method is suitable to select the largest value between 
positive slope up-crossings of the threshold and to eliminate the majority of smaller 
amplitude extremes that are less significant for the current purpose. The choice of 
the threshold is pretty important to investigate the statistics of the response and 
several methods are available [90] but since it is not the main point here, it is just 
used the threshold recommended by IEC61400-1, annex F, which is the mean value 
of the original series plus 1.4 the standard deviation. 

Figure 6.4 compares the exceeding probability distributions for both the model 
taking into account the wave impacts and the model considering only Morison’s 
contributions (solid red line and solid blue line, respectively). As expected, the upper 
tail of the impact model distribution keeps always above the distribution without 
impacts. If the distribution of the structural ultimate strength were also plotted in 
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Figure 6.2: Comparison of tower base overturning moment M,t due to EWM plus a 
JONSWAP irregular sea with and without the impulsive contributions due to breaking 
waves. 
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Figure 6.3: Comparison of tower base overturning moment M,; due to EWM plus a 
JONSWAP irregular sea with and without the impulsive contributions due to breaking 
waves. 


the same figure, it would result clear that the distribution without impacts leads 
to an overestimation of the structural capabilities giving rise to unsafe conclusions 
when assessing the probability of failure. 

The empirical exceeding probability distributions of peaks shown in figure 6.3 
is plotted in figure 6.4. It is fitted with a Generalized Extreme Value Distribution 
(GEV) as presented in figure 6.5. 

For the sake of simplicity we can imagine that the ultimate moment strength 
Mý: the tower can supply is a deterministic value. So we can draw a vertical line 
intercepting the two distributions plotted in figures 6.4 or 6.5 in two points repre- 
senting the probabilities that such an ultimate value is not exceeded for the model 
with impacts and for the model without impacts, respectively. These two values 
prove that the model without impacts leads to relay on a certain level of structural 
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Figure 6.4: Comparison of the exceeding probability distributions of the tower base 
bending moment M,; with and without considering the impulsive forces due to over- 


turning plunging breakers. 
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Figure 6.5: Comparison of the fitted GEV exceeding probability distributions of the 
tower base bending moment M,; with and without considering the impulsive forces 


due to overturning plunging breakers. 


safety which is actually wrong. In other words, the structural safety is overestimated 
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because the real probability that My, is not exceeded is given by the distribution 
accounting for the wave impacts. On the other hand, in the opposite situation, let 
us say in a design phase, one could fix a certain structural reliability which in this 
case can be represented by a horizontal line drawn on figure 6.4. The model not 
accounting for the impacts would lead to demand a lower ultimate bending strength 
than the more realistic one described by the distribution accounting for the impacts. 
Hence, we have that at 


Design stage: Given the Prai the non-impact wave model will always demand a 
lower design load; 


Verification stage Given the ultimate structural strength M},, the non-impact 
wave model will always underestimate the Pyai of the System. 


The above remarks are extremely relevant from a qualitative point of view and 
they pave the way for the development of a systematic methodology in assessing 
the ultimate loads for offshore wind turbines by employing the proposed model. 
However, in this context, since the probability distributions have been obtained 
only from one simulation, a reliable quantification of the error in the assessment of 
the probability of failure due to the inaccuracy of the model not accounting for the 
impacts cannot be provided. 


6.2 Summary and conclusions 


At the very end of this work it remains to draw a balance of what was planned 
to achieve within this research and what has been really attained. 

We started in Chapter 1 with a general overview about the European and world 
wind energy scenario. In this framework the most pressing priorities established by 
the European Commission to meet the challenging targets in terms of offshore wind 
energy have been recalled in order to provide some undiscussed motivations of gen- 
eral interest. In addition to them, further motivations have been identified by discov- 
ering some important lacks in modeling offshore wind turbines. In particular, while 
for fatigue design of offshore energy converters the current approach is satisfactory, 
when dealing with extreme wind and waves conditions very dangerous consequence 
in evaluating the structural safety may be induced by adopting standard methods 
alone. Departing from this background we theorized that an improved numerical 
model capable of capturing more realistically the response of a system exposed to 
an extreme environment could lead to a more accurate structural risk assessment. 
To this aim, after devoting some attention to the basic concepts concerning the 
rotor aerodynamics in Chapter 3, a fully nonlinear water waves simulator has been 
developed in Chapter 4. The code, written entirely within this research work, im- 
plements a higher-order Boundary Element Methods (BEM) to discretize in space 
Laplace’s equation governing the gravity waves propagation. Then, integrating in 
time the kinematic and dynamic boundary conditions - following indeed the so called 
Mixed Eulerian-Lagrangian scheme -, it is able to simulate fully nonlinear waves 
generated by several types of initial and boundary conditions. A numerical wave 
tank, equipped with an absorbing beach, has also been simulated and comparisons 
with experimental measurements, along with analytical results (when available), 
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validated the numerical tool. The Boundary Element Method-based code permits 
to reproduce plunging breaking waves up to the time the tongue re-enters into the 
free surface. 

The capabilities of the software developed in Chapter 4 have been used to set 
up an impact wave model which in turn has been embedded into a stochastic en- 
vironment. Indeed, the first part of Chapter 5 is devoted to the implementation of 
a logic that, once an extreme random sea is generated, whenever breaking waves 
occur, the fully nonlinear code is called to perform a refined analysis over a pace- 
time sub-domain defined as a neighborhood centered in the wind turbine location 
at the time of the expected breaking wave. Next, in the second part of Chapter 5, 
forces stemming from overturning breaking waves are passed to the hydro-aero- 
elastic simulator which permits to simulate the whole offshore wind turbine system 
when exposed to severe conditions. Results of such coupled model have then been 
thoroughly commented throughout the chapter. More synthetically, in this thesis: 


e A new numerical procedure to simulate extreme response of Offshore Wind 
Turbines has been developed: 


— The BEM-based code reproduces with high accuracy the overturning 
plunging breakers; 


— An analytical impact model permits to compute the impulsive forces 
subsequently passed to the aerolelastic solver; 


— The deterministic simulator has been successfully embedded into a stochas- 
tic environment; 


e Peaks in the structural response due to slapping waves can be up to three 
times higher than peaks induced by the standard linear wave approach. 


Finally, in this conclusive chapter we came back to what theorized at the begin- 
ning of the thesis, and in particular we tried to face those issues raised in Chapter 2. 
Here it has been shown that the new numerical model meets all the expectations. 
In fact: 


e The model developed significantly contributes to turn the simulation of off- 
shore wind turbines exposed to random wave impacts into a mature stage; 


e It represents a more reliable tool in predicting ultimate loading conditions; 


e It poses a crucial issue about the accuracy in evaluating the structural safety 
when non-impact models are employed because dangerous impact waves would 
be completely missed; 


e The model developed aims at being a design tool as it increases the model 
accuracy without penalizing the computational cost normally required; 


e An accurate extreme value analysis is necessary for a full quantification of the 
Structural Vulnerability. 


AppendixA 
Linear wave formulas 


Wind generated waves can be grouped in two main categories: 


Deep water waves, also referred to as short waves. The water is considered deep 
if the water depth d is larger than half the wavelength L, so d > L/2. These 
relatively short waves do not “feel” the sea floor. 


Shallow water waves, also known as long waves. The water is considered to be 
shallow if the water depth d is less than 1/20 of the wavelength. Namely, 
d < L/20. In this case the sea bottom has a significant influence on the 
characteristics of these relatively short waves. 


Thus, in the case of shallow water conditions (finite water depth) the kinematic 
and dynamic quantities for a single-harmonic (regular) wave are given as follows: 


Velocity potential 
_ gacoshk (y+ d) 


WwW cosh kd 


o (p, t) in (ka — wt) (A.1) 


Dispersion relation 
w? = gk tanh (kd) (A.2) 


Wavelength-wave period relation 


L=T /gd (A.3) 


Free surface elevation 


n (x,t) = acos (ka — wt) (A.4) 
Dynamic pressure 
cosh k (y + d) 
t)= | ka — wt A. 
pp (p,t) = —pgy + pga— cos (kx — wt) (A.5) 
Particle displacement along the x-axis 
___coshk(y+d) . 
¿E (p, t) = —a smhkd 22 (ka — wt) (A.6) 
Particle displacement along the y-axis 
sinh k (y + d) 
t) = a— ka — wt A. 
n (p,e) = YO cos (ka — wt) (A-7) 
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Velocity along the x-axis 


x ____ coshk (y+ d) 
v” (p,t) = aw kd 608 (ka — wt) 
Velocity along the y-axis 
J a sinh k (y +d) . E 
v” (p, t) ai n (ka — wt) 
Acceleration along the x-axis 
hk d 
a” (p,t) = qr Ok +d) sin (ka — wt) 


Acceleration along the y-axis 


2 sinh k (y + d) 
Oe sinh kd 


a” (p,t) = 


cos (kx — wt) 


(A.8) 


(A.9) 


(A.10) 


(A.11) 


Whereas, in deep water conditions (infinite water depth), the kinematic and 
dynamic quantities for a single-harmonic (regular) wave are given as follows: 


Velocity potential 
d (p, t) = e* sin (ka — wt) 


Dispersion relation 

Wavelength-wave period relation 
T 
On 


Free surface elevation 
n (x,t) = acos (ka — wt) 


Dynamic pressure 
pp (p,t) = —pgy + gae™” cos (ka — wt) 


Particle displacement along the x-axis 

E (p, t) = —ae*Y sin (ka — wt) 
Particle displacement along the y-axis 

n (p, t) = ae” cos (kx — wt) 
Velocity along the x-axis 


v” (p, t) = awe cos (kx — wt) 


(A.12) 


(A.13) 


(A.14) 


(A.15) 


(A.16) 


(A.17) 


(A.18) 


(A.19) 
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Velocity along the y-axis 
v” (p, t) = awe" sin (ka — wt) (A.20) 


Acceleration along the x-axis 


2 


a” (p,t) = awe" sin (ka — wt) (A.21) 


Acceleration along the y-axis 


a” (p, t) = —aw?e*Y cos (ka — wt) (A.22) 


where a = H/2 is the wave amplitude, k is the wave number, w is the circular 
frequency, d the water depth. The coordinates (x,y) of points p are referred to a 
Cartesian system having the y-axis upwardly oriented and the x axis positive in 
the wave propagation direction with the origin on the still water level. 


AppendixB 
Numerical dicretization of Laplace’s equation 


B.1  Green’s formulation 


Laplace’s equation (4.8) discussed in Chapter 4 together with the boundary and 
initial conditions represents a Boundary Value Problem (BVP) whose approximate 
solution is obtained by means of weighted residual methods [38]. The weighting 
function ¢* (p, pe) is the Green function that, for a 2D problem, reads as follows 


1 1 


®* (p, Pc) = Hi In R (B.1) 


where R = (a — xp.) + (Yp — Yp.)” is the distance between the point p and 
the collocation point pe. By means of equation (B.1) as well as the first and second 


Green’s identities, the BVP is turned into the following Boundary Integral Equation 
(BIE). 


cP) (pe) + fe WE f 6 (v.02) a(p) at = 0 (B.2) 


where 


a(p)=V@(p)-n (B.3) 
q” (p, pe) = V@* (p, pe) © (B.4) 


are respectively the flux (normal velocity component) and the normal derivative of 
the Green function. The coefficient c (pe) depends on the positions of the collocation 
point pe. In general it is given as 0/27, where J denotes the internal angle in radians 
of the corner at pe. When there is no corner, that is when V = 7, i.e. when pe lies 
on a smooth boundary, then c = 1/2; when the collocation point is internal to the 
fluid domain c = 1 and finally when pe is outside the domain c = 0 [38]. 

The boundary integral equation (B.2) is discretized into NE isoparametric qua- 
dratic elements in such a way nodal values belonging to j-th element are approxi- 
mated as follows 


3 

pD =F vx (8) of? (B.5) 
k=1 
3 

qP =Y ve (s) qe? (B.6) 
k=1 
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so that equation (B.2) becomes 


c (pe) È (Pe) Dai al Dg Dar 


Ti k=1 
E Di Sve 9A (B.7) 
Ti k=1 
where | | | 
pO = EPT) ns!) + (49 — yp.) ra” (B.8) 
i Fi [e — tp.) + (VO — Up.) | 
and 
; 2 5 
HO [0a mJ] a 


The shape functions adopted in equation (B.7) are as follows 


01(9)=356-1); 92()=(1-9(1+5); gs(s)=59(1+8) (B-10) 


Integrals in equation (B.7) can be transformed into integrals over the shape 
functions domain [-1, 1] as follows 


c (pe) $ (Pe HÈ Dal DO (s) TO (s)ds+ 


~ > NA > vr (8) A br (s) TO (s)ds = 0 (B.11) 


If the collocation point pe assumes in turn all NN boundary nodes, the above 
equation becomes 


NE 3 NE 3 
cidi + YO YO APP YI VP = i=1: NN (B.12) 
j=1 k=1 j=l k=1 


that is 


cipi + ADO? + DOS? + ALP GY? +... + AQP (NP) + AQP El” + ng gO”) + 
1 1 1 1 1 1 NE NE NE NE NE NE 
= PO? — go — gf OY? -—...- gli Pg) — gD” Of — GP of” =0 


where it has been set 
e Ra l 
R= | EO (r(e) T(s)as (B.13) 
=j 


tee T 8; (8) pr (8) TO (s)ds (B.14) 
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The number NN of boundary nodes is given by NN = 2NE+number of corners 
of the domain Q. 

To compute the velocity of particles internal to the fluid domain, the gradient 
of the velocity potential is obtained by differentiating equation (B.12). Thus, if the 
collocation point p, assumes in turn all NN;ng where the velocity is sought, the 
gradient of the velocity potential writes as follows 


NE 3 NE 3 . | 
P= ID Ga — YOVOAD dP -ASi NNint (B.15) 
j=1 na k=1 
NE 3 
cap sal ae yr: O G21: NNint (B.16) 
j=1 k= gel'k=1 
where 
. 1 
9, = / SO (5) pr (s) TO (ds (B.17) 
G) ta) 
Git =f Di? (8) pr (8) J (s)ds (B.18) 
1 
hive = g ) pr (8) J (s)ds (B.19) 
1 
ne al 2 (a ) ex (5) J (s)ds (B.20) 


For the sake of completeness, the gradient of the fundamental solution and its 
normal derivative for the i-th internal point are also given in the following 


sid 2 (2) = zi) 
bie ea (0 _ Li = (y) = TE (B.21) 
«gi l 2 (yP) — yi) 
iy Ch (O — n Ors ae (B.22) 
and 
( 0) i). + (09) la [(e = 21) nP + (y= y) nP] (2 — 2) 
*(7 4 
a 2r (CO + (yO yi)" ] 
i (B.23) 
P [( G) i) + (y vi) Jy- [(2 — a)n + (v= un n|(v— y) 
ned On (CO) ~ 2)? + (yO- 1] a 


The Jacobian 79) (s) referred 
differentiating the shape functions 


T9 (8) = NI 


(B.24) 


to the j-th boundary element is evaluated by 
as follows 
i 


J 


dzi) 
ds 


dy@) 


(B.25) 
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where the above derivatives are computed as follows 


de®) dog 
zO _ yo der) 


ds ads * (B-26) 
; 3 
= =>: Sek yl? (B.27) 
Finally, equation (B.12) can be rewritten in matrix form as follows 
CiPi + Pigha = Gigdq (B.28) 
for which 
Higha = Giga (B.29) 


where Hig = Cidig + Hig andei=1: NN and q=1:3NE. 
In the same way, equations (B.15) and (B.16) can be written in matrix form as 
follows 


vi =Gig,x4q = ig. Pq (B.30) 
uy =Gig.ydq — Hig Oa (B.31) 


where i = 1: NNint and q = 1:3NE. 

As soon as the primary unknowns of the system are computed (see the next two 
sections) the gradient of the velocity potential expressed by equations (B.30) and 
(B.31) can be readily computed. 


B.1.1 Assembling 


The matrix system (B.29), as well as (B.30) and (B.31), are of course not square 
because each boundary element has been treated as a single part without being 
connected with its two adjacent elements. An assembling algorithm is therefore 
required. The connectivity matrix has the following form 


Element nl n2 n3 

1 ni) ny? ns? 

de sE: SÙ si 

5 2 j-1 ae ne? ng» 
onn = . (3) (3) (3) 

J Ti n2 n3 
j+1 nt) nit) nth) 
NE Mi no am) 


Once the connectivity matrix has been defined, algorithm 4 turns the system 
matrices H and G into the square form HH and GG, respectively. By means of the 
same algorithm, matrices H z, H y, G x and G y are transformed in the square forms 
denoted by HH z, HH y, GG and GG y. 
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Algorithm 4: Assembling 


input : H and G of size NN x 3NE 
output: HH and GG of size NN x NN 


for i=1to NN do 
for j = 1 to NE do 
for k = 2 to 4 do 
index = Conn (j, k) 
q = 3j —34+(k-1) 
HH (i, index) = HH (i, index) + H (i,q) 
GG (i, index) = GG (i, index) + G (i,q) 


B.1.2 Reordering the system and continuity conditions 


The system 
HH;;®; = GGijqj (B.32) 
needs to be rearranged according to the boundary condition type in order to get an 
algebraic system in the standard form 


Ais Xj = bi (B.33) 


where all the unknown are collected in X and all the known terms, among velocity 
potential and flux, are collected in b. 

Singularities at the upper left and upper right corners of the domain, where two 
different boundary conditions coexist, have been removed by the so called “double— 
node technique”. However, this expedient causes that for each of the two corners 
the matrix A has two identical row and consequently the system results singular. 
To avoid this the continuity condition is inserted according to [43] (see section 4.3 
of the paper). 

Finally, once the system unknowns are found, the gradient of the velocity po- 
tential at points internal to the fluid domain is computed as follows 


a 


vi =GGiq,y4j == AE gy; (B.35) 


a 


where 1 = 1: NNint and j= 1: NN. 
B.2 Numerical dervatives 


Once the velocity potential is known on the free surface, the tangential compo- 
nent of the particle velocity at each node is computed by the following numerical 
scheme 

vy =VO-t (B.36) 


It is straightforward to prove that the above equation can be rewritten as follows 


i (ddr | ög dy 1 dọ 1 
T= È ds dy si) Ta ds J(s) (B 
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where for the j-th element the velocity potential derivative with respect to the 
curvilinear coordinate s, as well as the derivatives involved in the Jacobian, are 
obtained by differentiating the shape functions as follows 


dọ _dgi (j) 

EA REA B.38 
ds ds di ( ) 
dz dy; G) 

pentyl N B. 
ds ds? (8:39) 
dy _ dei) 

Re IE A B.40 
ds ds Yi ( ) 


B.3 Gradient of the velocity field 
Let us start form the gradient of the velocity field 
vi = ( Dia Wry ) 
Vry  Vyy 
The curve Ty : R —> IR?, i.e. the free surface, is parameterized as follows 


sH T; (8) “i als) 
The velocity field is also parameterized as follows 


tee 2 f Oe (Ets Gs) st) 
sH U(s) -{ vY (x (s),y(s),t) 


The chain rule allows to write the following partial derivatives 


dv? dv der ðv” dy 


= + B.41 
Os dx ds dy ds ( ) 
ðv! Ov¥ dx dv dy 
= B.42 
Os dx ds dy ds ( l 
Next, the continuity equation assures that 
dv” dv 
=- B.43 
Ox Oy ( ) 
while the irrotational flow hypothesis allows to write 
ðv! ov? 
—— = — B.44 
Ox Oy ( ) 


and both of the above equations transform the system of equations (B.41) and 
(B.42) into 


ðv” Ov der ðv” dy 
ðs Oxds Ox ds 
dv dv dx | Ov? dy 
ðs Oxds dx ds 


(B.45) 


(B.46) 
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which can be solved for ee and dei to get 
dv” 1/0v da dv dy 
= B.4 
ðr D ( Os ds ðs =) CEI 
ðv” 1/0vdy dvda 
ðr D ( ðs ds = ds =) (B48) 
(B.49) 


where 


D(s) = J (s)? = (=). n eh (B.50) 


Finally the second term of the right-hand side of equation (4.21) becomes 


D T 

n = OF HVF UY + UF ve (B.51) 
Dv. 

pi” OF + UG 0° HUG yV” (B.52) 


where the components on and on have been found by using the continuity and 


irrotational equations. 
B.4 Tests of convergence 


In order to check the reliability of the BEM code it is important first to show the 
reliability of the solution of the steady Laplace’s equation. This is in fact the core 
of the numerical wave simulator because at each time step a stationary problem is 
solved. 

To this aim in the first two tests a very well known problem in solid mechanics 
is presented: the Saint Venant torsional problem. Indeed, given a cylindrical beam, 
i. e. a solid with one dimension dominating with respect to the others, and denoting 
its cross section by A, the whole state of stress induced by a torsional couple is 
known when either the stress function y or the warping function £ are known. 
Details on the torsion problem are available in [92]. 

The last case regards a a general potential problem defined on a square domain 
with mixed boundary condition of the type as in the case of water waves. 


B.4.1 Torsion: Dirichlet’s problem 


The problem can be formulated in terms of the stress function w as follows 


V?24 (21, £2) =0 VpEA 
ili VpedA (B-38) 


When the cross section A has an elliptical shape, the analytical solution writes 
as follows sea PERI 
afb laf — b 
= H B.54 
gui 9) oe 


where a and b are the major and minor axes of the ellipse. 
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NE | 2 4 8 16 32 64 


err L2 | 74.9129 0.7892 0.0017 0.0000 0.0000 0.0000 


errLo | 8.6552 0.8884 0.0408 0.0028 0.0001 0.0000 
Table B.1: errL3 and errL of stress function for different boundary mesh. 
Given such a domain, for four boundary discretizations more and more refined 


we compare the numerical and analytical solutions by using the norm Lə of the 
error computed as 


n 


errLo = y (pára. — gNum.)? (B.55) 


i=l 


where n is the number of points uniformly distributed across the domain A at 
which the two solutions are evaluated. See the red dots into the domains plotted 
in figure B.1. The figure shows different boundary discretizations form a number of 
boundary elements NE = 2 up to NE = 64. 

The error decrease quadratically and approaches zero with sixteen elements. See 
figure B.2. The rate of convergence is pretty satisfactory. 


B.4.2 Torsion: Neumann’s problem 


The excellent result shown in the previous section needs to be confirmed also in 
the case when mixed boundary conditions are applied (as the case of our numerical 
wave simulator) where discontinuities occur. As already mentioned, this problem is 
faced by introducing the double-node technique. To test such a circumstance, the 
domain is here restricted only to one quarter of the all ellipse. 


(B.56) 


V7 (x1, £2) =0 VpeAA 
Vo: n = TNn — CM VpedA 


In the case of Dirichlet’s problem no corners were involved, thus the double-node 
method was not necessary. On the contrary, for the one forth of ellipse, there are 
three corners. For those where Neumann-Dirichlet boundary conditions are assigned 
the double-node technique works quite well after implementing the the kinematic 
and continuity conditions proposed in [43]. But at the symmetry point, in the center 
of the ellipse, the Dirichlet-Dirichlet condition produces two identical equations. 
Hence the system is again singular. To avoid this, the two nodes have been shifted 
of a very small quantity inside the elements [38]. With this expedient results are 
sufficiently accurate. 

The domain discretization is shown in figure B.3. While the Root Mean Square 
Error convergence in shown in table B.2 and figure B.4. 


B.4.3 Potential problem on a square domain 


In this case a square domain Q whose boundaries Ty, r4, T f, Ti are straigth lines 
with lengths Ly = Li, Ly = Li, 1m is considered. The Neumann boundary 
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BEM solution of S. Venant Torsional Problem, NE = 4 
T T £ 


4 | of | 
ry Tı 
(a) NE=2 (b) NE =4 
BEM solution of S. Venant Torsional Problem, NE = 32 BEM solution of S. Venant Torsional Problem, NE = 64 


(c) NE = 32 


(d) NE = 64 


Figure B.1: Mesh refinement in the case of an elliptical cross section, domain for a 


Dirichlet’s torsional problem. 


0 10 20 


30 


NE 


50 60 70 


Figure B.2: Convergence of errLo in the case of Dirichlet’s torsion problem. 
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BEM solution of S. Venant Torsional Problem, NE, =32, NE. =1, NE» =16 BEM solution of S. Venant Torsional Problem, NE =32, NE. =32, NE2 =16 


(a) NE, NEe, NE3 = 32, 1,16 (b) NE, NE», N Ez = 32,128, 16 


Figure B.3: Mesh refinement in the case of an elliptical cross section, domain for a 
Neumann’s torsional problem. 


0 20 40 60 80 100 120 140 


Figure B.4: Convergence of errLo in the case of Neumann’s torsion problem. 
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NEe | 1 2 4 8 16 32 64 128 


err L2 | 2.5598 0.2588 0.0183 0.0029 0.0008 0.0003 0.0002 0.0002 


errL2 | 1.5999 0.5087 0.1353 0.0539 0.0278 0.0177 0.0139 0.0125 


Table B.2: errL3 and errL2 of torsion function for different boundary mesh. 


Li ¢ 
-0.3 -0.3 
-0.4 -0.4 
D 0.5) € ¢ > -0.5 
-0.6} -0.6 
-0.7 -0.7 

e e 


-1p € . Li . € | iL 


0 0.1 02 03 04 05 06 07 08 09 1 0 0.1 02 03 04 05 06 07 08 09 1 
T z 


(a) Number of elements per side NE =2 (b) Number of elements per side NE = 16 


Figure B.5: Two discretizations for the solution of the potential problem defined on 
a square domain. 


conditions are prescribed as follows 


acos (a (y + d)) Vp ELl, 
200) = —ae~*cos(a(yt+d)) Vper,, 
n 0 Vp € Ty 


while on the remaining boundary the Dirichlet condition is ® (p) = 0 Vp € Ty. a is 
a constant equals to 37/2 and d = Li, = Li, is the distance between the upper and 
lower boundaries (in the next cases it will be the water depth). The exact solution 
of this problem is ® (p) = e~** cos (a (y + d)). 

The vector X collecting all the unknowns of the problem concatenates 2NE, + 1 
elements that are the potential function over Te, 2NE;, + 1 elements that are the 
potential function over [;,, 2NE;+1 elements that are the fluxes over T ș and finally 
2NE;, + 1 elements that are the potential function over T;. 

To investigate the error convergence of the steady second order boundary ele- 
ment solver, the global error errL2 has been evaluated as function of the number 
of elements per boundary. All the four boundaries have been discretized with the 
same number of elements. See figure B.5. 

Figure B.6 shows the convergence of the error computed as 


errL9 = (B.57) 
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Error vs. NE 


n i i 
10° 10° 10° 10 


Number of elements per side NE 


Figure B.6: Error convergence of the error. 


where X and X, are the vectors collecting the numerical and the exact analytical 
solution, respectively, only on the boundary. 

The global error decreases quadratically as the number of element per boundary 
NE increases. This rate of convergence is justified by the order of the boundary 
elements discussed above. When the number of elements per side is quite low, when 
it doubles from 2 to 4, the convergence is even faster then quadratic. 
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